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Abstract 

Accurate  change  detection  (CD)  results  in  urban  environments  is  of  interest  to  a 
diverse  set  of  applications  including  military  surveillance,  environmental  monitoring, 
and  urban  development.  Although  CD  approaches  with  hyperspectral  imagery  exist, 
there  is  no  comprehensive  framework  that  describes  existing  approaches  and  provides 
a  foundation  to  generate  new  ones.  This  work  presents  a  hyperspectral  CD  (HSCD) 
framework.  The  framework  uncovers  the  need  for  HSCD  methods  that  resolve  false 
change  caused  by  image  parallax.  Image  parallax  is  the  apparent  motion  of  stationary 
objects  due  to  differing  viewing  geometries. 

A  Generalized  Likelihood  Ratio  Test  (GLRT)  statistic  for  HSCD  is  developed 
that  accommodates  unknown  mis-registration  between  imagery  described  by  a  prior 
probability  density  function  for  spatial  mis-registration.  Using  a  normal  prior  dis¬ 
tribution  for  the  mis-registration  leads  to  a  fourth  order  polynomial  test  statistic 
that  can  be  numerically  minimized  over  the  unknown  mis-registration  parameters.  A 
more  computationally  efficient  closed-form  solution  is  developed  based  on  a  quadratic 
approximation  derived  through  a  second  order  Taylor  series  expansion  and  provides 
comparable  results  to  the  numerical  minimization  for  the  investigated  test  cases  while 
running  30  times  faster.  The  results  indicate  an  order  of  magnitude  reduction  in  false 
alarms  at  the  same  detection  rate  for  synthetically  mis-registered  test  data,  especially 
in  image  regions  containing  edges  and  hne  spatial  features.  Investigation  of  the  model 
parameters  is  also  assessed,  and  the  potential  of  the  derived  method  to  incorporate 
more  complex  signal  proccessing  functions  is  demonstrated  by  the  incorporation  of  a 
parallax  error  mitigation  component.  The  GLRT  statistic  is  extended  to  account  for 
parallax  errors  and  the  visual  results  demonstrate  the  reduction  of  false  alarms;  also 
evident  in  test  statistic  histograms  and  signal-to-clutter  ratio  (SCR)  computations. 
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REMOVING  PARALLAX-INDUCED  FALSE  CHANGES  IN  CHANGE 


DETECTION 


I.  Introduction 


1.1  Overview 

Change  detection  (CD)  is  the  problem  of  discriminating  significant  changes  from 
insignificant  changes  between  two  images  [6].  Traditional  image-based  CD  methods 
nse  red-green-blne  color  imagery.  Dne  to  the  spectral  limitations  with  these  imaging 
technologies,  the  remote  sensing  commnnity  often  exploits  hyperspectral  (HS)  im¬ 
agery  for  the  CD  task.  A  literature  review  uncovers  relatively  few  works  that  develop 
novel  CD  approaches  for  use  on  HS  imagery.  Certainly  there  exist  benefits  of  per¬ 
forming  HSCD  over  regular  image-based  CD,  such  as  detecting  imperceptible  targets 
in  a  complex  background.  However,  challenges  both  in  integrating  the  spectra  into 
the  change  detector  and  working  with  imagery  derived  by  non-stationary  platforms 
makes  HSCD  particularly  difficult. 

Accurate  CD  in  urban  environments  is  of  interest  to  a  diverse  range  of  communities 
and  a  diverse  range  of  applications  to  include:  military  surveillance,  environmental 
monitoring,  and  urban  development.  Manually  processing  of  data  for  CD  is  daunting; 
especially  in  the  case  of  HS  data,  due  to  its  hundreds  of  spectral  channels.  Therefore, 
there  is  a  strong  need  for  methodologies  that  enable  automated  detection  of  changed 
or  interesting  events  in  aerial  imagery  to  alert  an  image  analyst  or  operator.  Such  a 
methodology  could  provide  further  support  by  identifying  relevant  regions  of  change. 
The  incentive  for  performing  HSCD  is  to  take  advantage  of  all  the  possible  change 


1 


information  available  in  the  data  providing  a  nniqne  method  for  identifying  snbtle 
changes  in  a  complex  environment. 

1.2  Research  question 

A  basic  CD  algorithm  generates  a  change  mask  that  identifies  change  between 
pixels  from  a  set  of  images  acquired  at  different  times.  The  reliability  of  the  change 
mask  is  influenced  by  factors  such  as  camera  motion,  sensor  noise,  illumination  vari¬ 
ation,  non-uniform  attenuation,  or  atmospheric  absorption.  Radke  et  al.  [B]  illustrate 
that  detecting  changes  and  determining  their  signihcance  is  complex  and  challenging. 
Two  common  problems  with  the  application  of  CD  algorithms  for  airborne  data  are 
image  registration  and  parallax  error  [7]  between  the  reference  and  test  scenes  caused 
by  different  viewing  geometries  between  the  images.  Registration  errors  between  the 
image  pair  causes  pixel  mis-alignment  (mis-registration)  and  impacts  pixel-based  CD. 
Parallax  error,  the  apparent  change  in  relative  positions  of  stationary  objects  due  to 
different  viewing  geometries,  is  prevalent  with  taller  objects  in-scene  and  often  causes 
false  change  [H].  Image  mis-alignment  further  impacts  the  ability  to  determine  if  im¬ 
age  parallax  is  the  cause  of  false  detections  while  image  parallax  makes  the  alignment 
of  previously  present  objects  impossible,  also  resulting  in  false  detections. 

The  signihcant  impact  of  image  mis-registration  and  image  parallax  on  HSCD 
leads  to  the  following  research  question:  Can  one  develop  a  HSCD  method  that  si¬ 
multaneously  removes  false  changes  caused  by  image  mis-registration  and  parallax? 

1.3  Methodology  preface 

The  theory  and  phenomenology  of  hyperspectral  processes  that  support  change 
detection  is  introduced.  Specihcally,  derivations  of  common  change  detectors  for  hy¬ 
perspectral  data,  chronochrome  (CC)  and  covariance  equalization  (CE)  are  provided 
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as  current  state-of-the-art.  A  new  method  to  account  for  image  mis- registration  and 
parallax-induced  false  detections  in  HSCD  is  presented.  The  general  form  of  the  bilin¬ 
ear  interpolation  equation  is  used  to  account  for  estimates  of  image  mis-registration 
for  the  sub-pixel  case,  where  a  generalized  likelihood  ratio  test  (GLRT)  is  constructed 
to  provide  a  change  test  statistic  for  change  at  each  pixel  in  the  temporal  hyperspec- 
tral  image  pair.  This  approach  is  then  modified  to  incorporate  an  estimate  of  the 
parallax  direction  in  order  to  mitigate  false  alarms  due  to  parallax  affects  between 
the  same  two  hyperspectral  image  pairs. 

1.4  Results  preface 

After  outlining  the  practical  implementation  and  data  sets  used,  applications  using 
hyperspectral  data  are  used  to  detail  the  GLRT-based  algorithms  presented  and  are 
compared  with  current  anomalous  change  detectors.  The  results  of  the  implementa¬ 
tions  are  discussed  and  performance  of  the  GLRT-based  algorithms  is  evaluated  using 
receiver  operating  characteristic  (ROG)  curves  when  applicable.  Gomparisons  among 
curves  are  evaluated  in  terms  of  area  under  the  ROG  curve  (AUG)  estimated  on  the 
basis  of  available  ground  truth.  The  AUG  for  the  mis-registration  GLRT-based  al¬ 
gorithm  (AUGglrtcc  =  0.9981  and  AUGglrtcb  =  0.9954)  is  higher  than  the  current 
state-of-the-art  baseline  algorithms  (AUGcc  =  0.9972  and  AUGge  =  0.9909)  which 
demonstrate  the  effectiveness  of  the  methods  established  in  this  research.  In  relation 
to  a  similar  mis-registration  approach.  Local  Go- Registration  Adjustment  (LGRA), 
the  GLRT-based  algorithm  provides  comparable  performance.  Visually,  test  statistics 
reveal  that  parallax  mitigation  reduces  false  alarms. 
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1.5  Organization  of  dissertation 


Chapter  presents  the  common  processes  used  in  CD  and  computer  vision  per¬ 
tinent  to  this  research.  A  process  framework  for  HSCD  is  introduced  in  order  to 
frame  the  theoretical  developments  of  this  work.  The  new  theoretical  foundation  for 
addressing  false  detections  due  to  image  mis-registration  is  described.  Furthermore, 
Chapter  [ITT]  builds  on  this  foundation  to  mitigate  false  detections  due  to  image  paral¬ 
lax.  Chapter  |IV|  provides  detailed  simulation  results  demonstrating  the  effectiveness 
of  the  mitigation  approaches.  Finally,  Chapter  |V|  summarizes  the  research  presented 
before  highlighting  the  contributions  and  efforts  for  potential  future  research. 
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II.  Background 


2.1  Introduction 

Detection  and  characterization  of  geospatially  coincident  and  temporally  variant 
imagery  has  been  of  interest  to  researchers  for  many  years  due  in  large  part  to  the 
number  of  applications  across  diverse  disciplines.  Automated  image  change  detection 
(CD)  is  the  process  of  autonomously  analyzing  spatially  coincident  image  pairs,  and 
identifying  regions  that  have  undergone  signihcant  spatial  or  spectral  variation.  CD 
is  an  important  process  in  monitoring  and  managing  military  applications  because  it 
provides  quantitative  analysis  of  the  spatial  distribution  of  the  population  of  interest. 

Macleod  and  Congation  [9j  list  four  aspects  of  CD  which  are  important  when 
performing  a  monitoring  task:  (1)  detecting  the  changes  that  have  occurred;  (2) 
identifying  the  nature  of  the  change;  (3)  measuring  the  spatial  area  extent  of  the 
change;  and  (4)  assessing  the  spatial  pattern  of  the  change.  CD  results  are  further 
processed  by  some  means  of  characterization  to  filter  the  number  of  detections  to  a 
more  meaningful  subset  based  upon  some  application  specific  requirement  or  a  priori 
information. 

The  basis  of  using  remote  sensing  data  for  CD  is  that  scene  changes  result  in 
radiance  value  differences  that  can  be  sensed  remotely.  Imaging  spectroscopy,  the 
process  of  simultaneously  acquiring  data  at  multiple  distinct  spectral  wavelengths,  is 
a  passive  remote  sensing  technology  that  is  of  specific  interest  to  the  study  of  ad¬ 
vancing  CD  capabilities.  Spectral  sensor-related  constraints  that  affect  the  quality  of 
the  imagery  and  bound  the  options  for  processing  image  data  include:  bandwidth, 
sampling  interval  (the  spacing  between  two  consecutive  wavelengths),  and  data  type 
(he.,  the  number  of  bits  used  to  digitally  represent  the  magnitude  of  the  sensor  re¬ 
sponse).  These,  and  other  constraints,  are  typically  referred  to  as  sensor  operating 
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conditions  (OCs)  [51 HH]-  Two  commonly  referenced  classes  of  spectral  imaging  spec¬ 
trometers  are  hyperspectral  (HS)  and  multispectral  (MS).  The  distinction  between 
these  two  classes  of  imagers  are  typically  discriminated  by  the  number  of  spectral 
bands  captured  as  well  as  the  channel  bandwidth;  an  acceptable  benchmark  for  HS 
sensors  is  the  collection  of  hundreds  of  contiguous  spectral  bands  while  MS  sensors 
collect  less  than  20  generally  non- contiguous  spectral  bands  (Fig.  [^.  An  additional 
discriminator  for  an  HS  versus  MS  imaging  system  is  the  full  width  at  half  maximum 
(FWHM).  FWHM  refers  to  the  detector  response  derived  from  exposure  to  a  source 
and  is  assumed  to  be  Gaussian.  The  FWHM  contributes  to  the  spectral  channel 
bandwidth,  typically  reported  in  microns  or  nanometers  (Fig.  i  mi- 

HS-based  CD  has  many  advantages  over  MS-based  CD  in  detecting  and  discrimi¬ 
nating  surface  materials  because  the  former  provides  a  continuous  spectrum  across  a 
wide  range  of  wavelengths  where  the  later  does  not.  Typical  gray-scale  imagery  (an 
image  in  which  the  value  of  each  pixel  is  a  single  intensity)  has  a  limited  capability 
for  detecting  spectral  change,  but  generally  provides  greater  spatial  resolution.  These 
gray-scale  images  often  fail  to  provide  sufficient  information  to  discriminate  low  con¬ 
trast  targets  that  might  be  employing  concealment  techniques.  Hyperspectral  change 
detection  (HSCD)  provides  the  added  sensitivity  to  reveal  camouflage,  and  in  general, 
camouflage,  concealment,  and  deception  (CC&D)  targets.  Here,  HSCD  is  addressed 
by  means  of  a  process  framework  developed  from  a  thorough  study  of  the  CD  liter¬ 
ature.  This  CD  framework  focuses  on  passive  remote  sensors,  in  particular  imaging 
spectrometers  that  only  measure  radiation  emitted  or  reflected  into  the  sensor  by  the 
object  under  scrutiny  and/or  its  surroundings. 


6 


400.0  700,0  10000  1300.0  1600.0  1900.0  2200.0  2500,0 

Wavelength  (nm) 


Figure  1.  A  comparison  of  the  NASA’s  Airborne  Visible/Infrared  Imaging  Spectrome¬ 
ter  (AVIRIS)  224  spectral  channel  sensor  for  a  measured  spectra  from  400  to  2500  nm 
at  10  nm  intervals  across  the  solar  reflected  spectrum  and  6  of  the  multispectral  bands 
measured  by  the  LANDSAT  Thematic  Mapper  [1]. 


Figure  2.  Illustration  of  full  width  at  half  marximum  (FWHM)  used  to  indicate  spectral 
resolution. 


2.2  Hyperspectral  change  detection  literatnre  review 

There  are  several  methods  for  addressing  HSCD  in  an  end-to-end  system;  these 
methods  are  organized  into  groups  for  examination:  pre-processing,  change  detectors, 
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thresholding,  post-processing,  and  analysis/assessment.  The  major  topics  are  exam¬ 
ined  to  provide  relevant  background  for  the  proposed  framework.  Although  there  is 
a  multitude  of  approaches  in  the  literature  as  to  the  organization  of  CD,  they  con¬ 
verge  on  the  following  recommendations:  CD  should  involve  data  acquired  by  the 
same  sensor,  with  the  same  spatial  resolution,  same  viewing  geometry,  same  spectral 
bands,  same  radiometric  resolution,  and  acquired  at  the  same  time  of  day  [7].  It  is 
also  desirable  to  perform  the  change  analysis  at  the  same  collection  time  to  limit  the 
differences  associated  with  sun  angle  and  seasonal  variations. 

The  simplest  CD  framework  separates  signihcant  and  non-signihcant  changes  that 
generate  a  correspondence  image  from  an  image  pair  (Fig.  [^.  Images  are  typically 
referred  to  as  the  reference  image,  the  test  image,  and  the  change  response.  In  the 
CD  routine,  two  corresponding  pixels  belonging  to  the  same  location  in  an  image  pair 
are  determined  on  the  basis  of  a  quantitative  measure.  Two  temporal  image  frames 
(or  a  pair  of  temporal  images)  are  the  minimum  requirement  to  identify  change. 


■=> 

Image  3 


Image  2 

Figure  3.  Simple  CD  framework  using  a  minimal  processing  schema.  Here  Image  1  is 
referred  to  as  the  reference  Image,  Image  2  the  test  image,  and  image  3  the  change 
response. 


Although  CD  can  be  performed  in  the  simplest  manner  as  depicted  in  Fig. 


there  is  evidence  that  an  extended  form  would  improve  HSCD  performance.  HSCD 
results  depend  on  many  factors,  and  this  section  details  these  factors  into  a  com¬ 
prehensive  framework.  Jensen  na  introduces  a  CD  procedure  to  monitor  land  CD 
using  remotely  sensed  data  that  consists  of:  (1)  state  the  CD  problem,  (2)  consid¬ 
erations  of  signihcance  when  performing  CD,  (3)  image  processing  of  remote  sensor 
data  to  extract  change  information,  (4)  quality  assurance  and  control  program,  and 
(5)  distribute  results.  Lu  et  al.  [13]  describe  three  major  steps  involved  when  imple¬ 
menting  a  CD  system:  (1)  image  pre-processing,  (2)  selection  of  suitable  techniques, 
and  (3)  accuracy  assessment.  Similarly,  Jianya  et  al.  [Tl|  discuss  the  problem  of  CD 
and  focus  on  the  same  three  steps  as  in  [T3|,  but  with  the  terms:  (1)  pre-processing, 
(2)  CD  algorithms,  and  (3)  accuracy  assessment.  These  simplistic  descriptions  of  the 
CD  processing  chain  lack  robustness  and  do  not  capture  many  of  the  challenges  in 
image-based  CD.  These  frameworks  further  do  not  take  advantage  of  HS  images. 

The  literature  is  in  common  agreement  that  a  CD  algorithm  requires  additional 
processes  beyond  the  simple  frameworks  described  previously.  However,  there  is  no 
consensus  as  to  a  single  universally  applicable  processing  chain.  The  choice  of  routines 
within  each  method  in  the  framework  must  be  navigated  based  on  the  application 
and  the  needs  of  the  end-user.  Once  the  areas  of  the  framework  are  determined,  an 
understanding  of  specihc  data  processing  requirements  is  easily  assessed. 

Relatively  few  works  support  HSCD  and  most  do  not  dehne  the  entire  processing 
chain.  Eismann  et  al.  naiiE]  propose  the  use  of  an  affine  transformation  to  eliminate 
errors  associated  by  multiple  algorithms  during  the  pre-processing  stage.  An  affine 
transformation  is  a  linear  transformation  and  translation  between  two  domains  -  for 
CD  this  is  the  reference  domain  to  the  test  domain.  Eismann  et  al.  [13  HB]  assume  the 
same  sensor  is  utilized,  images  are  acquired  under  similar  environmental  conditions, 
and  that  image  matching  or  cropping  has  already  occurred.  Ortiz-Rivera  et  al.  im 


9 


also  present  general  steps  for  HSCD  to  include  pre-processing  for  the  estimation 
of  the  change  mask,  which  indicates  the  pixels  that  reflect  a  change.  The  steps 
are  application  specihc  to  an  extent  and  should  be  further  developed  to  include  all 
pre-processing  possibilities  such  as  feature  extraction  and  dimensionality  reduction. 
Ultimately  the  environmental  setting  plays  an  important  role  in  the  ability  to  identify 
change;  as  the  setting  becomes  more  complex  (he.,  heavy  vegetation  surrounding  a 
city)  the  accuracies  of  CD  techniques  suffer.  The  literature  often  focuses  on  a  single 
physical  environment  limiting  the  impact  of  their  analysis  in  evaluating  different  CD 
methods  across  a  variety  of  background  conditions  PEHIE].  Incorporating  post¬ 
processing  into  the  CD  processing  chain  may  enable  a  more  diverse  study  by  level 
setting  results  due  to  various  background  conditions. 

2.3  Organization  of  CD  methods 

Table  provides  a  basis  for  discussing  a  unihed  CD  processing  chain.  It  lists  hve 
foundational  methods:  pre-processing,  change  detector,  post  processing,  thresholding, 
and  analysis/assessment  and  summarizes  typical  routines  within  each  overarching 
method  with  applicable  supporting  citations. 
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Table  1.  A  table  of  CD  methods  with  supporting  routines  and  literature  references. 


Change  Detection  Methods 

Pre- 

Change 

Post- 

Threshold 

Post- 

Analysis 

processing 

detector 

processing 

processing 

assessment 

•Geometric 

•Local- 

•HSV 

•Empirically 

•Ealse 

•Visual 

registration 

based 

Trans- 

•Theoretical 

alarm  mit- 

•Quantitative 

•Radiometric 

methods 

formation 

igation 

•ROC  curve 

correction 

•Pixel- 

•Contextual 

•Confusion 

(or  normal- 

based 

information 

matrix 

ization) 

methods 

Radke  [6] , 

Rosin  [23] , 

Radke  [6] , 

Jensen  [12], 

Knud- 

Sahoo  [25]. 

Radke  [6] , 

Rosin  [23] , 

Eis- 

Lu  [13], 

sen  [22]. 

Sezgin  [26] 

Lee  [27] 

Sezgin  [26] 

mann  [15] , 

Marce- 

Zhang  [23| 

Rosen- 

naro  |21j 

feld  [20] 

The  pre-processing  stage  adds  distinct  benefits  to  the  CD  process.  Precise  reg¬ 
istration  of  multi-temporal  imagery  is  a  critical  prerequisite  of  accurate  pixel-level 
CD.  The  change  detector  is  the  heart  of  the  change  detection  method  and  a  vari¬ 
ety  of  change  detectors  are  available.  A  threshold  scheme  is  typically  employed  to 
distinguish  a  separation  between  signihcant  and  insignihcant  changes.  The  selec¬ 
tion  of  an  appropriate  threshold  value  for  discriminating  unchanged  pixels  from  the 
changed  pixels  is  also  a  key  factor  in  most  CD  methods  |2H].  A  binary  change  mask 
is  created  by  classifying  those  pixels  with  a  value  above  the  threshold  as  being  a 
valid  change.  Unfortunately,  there  are  a  number  of  factors  that  make  obtaining  high 
accuracy  change  masks  a  difficult  task.  Two  common  examples  include  changing  il¬ 
lumination  conditions  and  image  pair  mis- alignment.  Post-processing  those  changes 
can  improve  the  change  mask  result. 

Post-processing,  though  rarely  mentioned  in  the  HSCD  literature,  is  used  to  refine 
CD  results  and  can  be  viewed  as  a  false-alarm  mitigation  method.  Rudimentary 
approaches  post-process  the  change  mask  with  standard  binary  image  processing 
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operations  [B] .  Other  approaches  build  binary  classifiers  based  on  training  data  from 
the  change  response  to  produce  a  more  accurate  change  result  [29] .  Table  includes 
post-processing  as  a  capability  before  or  after  the  thresholding  stage. 

Lu  et  al.  [T3|  provide  factors  that  influence  the  accuracy  of  CD  results,  including: 
(1)  precise  geometric  registration  between  multi-temporal  images,  (2)  calibration  or 
normalization  between  multi-temporal  images,  (3)  availability  of  quality  ground  truth, 
(4)  the  complexity  of  the  landscape  and  environment  of  the  study  area,  (5)  CD 
methods  or  algorithms  used,  (6)  classification  and  CD  schemes,  (7)  analyst’s  skills 
and  experience,  (8)  knowledge  and  familiarity  of  the  study  area,  and  (9)  time  and  cost 
restrictions.  The  performance  of  a  CD  algorithm,  including  the  necessary  pre-  and 
post-  processes,  can  be  evaluated  visually  and  quantitatively  depending  on  the  goal. 
Subcomponents  of  the  CD  framework  can  also  be  evaluated  individually  [6] .  Assessing 
the  accuracy  of  CD  products  is  an  important  step  in  assessing  the  validity  of  the 
method.  In  assessing  changes  based  on  remotely  sensed  data,  the  major  impediment  is 
that  the  estimate  values  are  difficult  to  compute  due  to  the  complexity  of  the  processes 
involved  and  more  often  the  truth  data  is  not  available  for  computing  accuracy.  When 
the  inherent  limitations  are  appropriately  dealt  with,  pre-processing  is  adequately 
incorporated,  and  appropriate  CD  algorithms  are  selected,  the  assessment  will  reveal 
adequate  results. 

Analysis  of  the  literature  provides  ample  evidence  to  support  the  conclusion  that 
HS  data  can  be  effectively  used  to  detect  and  monitor  changes.  However  there  is 
an  agreement  with  the  observation  of  [HU]  that  one  of  the  challenges  confronting  the 
community  is  to  develop  an  improved  understanding  of  the  CD  process  on  which  to 
build  an  understanding  of  how  to  match  applications  and  methods.  Our  CD  frame¬ 
work  provides  an  improved  understanding  based  on  careful  analysis  of  the  literature 
and  moves  the  community  closer  to  matching  applications  and  methods  for  improved 
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CD  results. 


2.4  Change  detection  framework 

The  framework  presented  in  this  dissertation  document  provides  a  foundation  for 
comparing  HSCD  systems  |2].  The  framework  provides  categories  of  methods  that 
must  be  acknowledged  when  analyzing  HSCD  and  is  presented  in  Fig.  It  begins 
with  two  temporal  HS  image  cubes,  a  reference  image  and  a  test  image.  The  first 
step  is  to  pre-process  both  image  cubes  to  suppress  undesired  distortions  and  enhance 
image  features  important  for  the  rest  of  the  processing  change.  The  pre-processed 
image  cubes  are  supplied  to  the  CD  step  which  then  outputs  a  CD  response  map 
with  the  raw  detector  results.  The  CD  response  map  can  either  be  post-processed 
or  proceed  to  a  threshold  step  in  order  to  identify  significant  changes  and  represent 
them  as  a  change  mask.  The  change  mask  can  then  be  enhanced  in  an  additional 
post-processing  step.  The  hnal  change  mask  is  then  used  to  analyze  CD  performance. 


Figure  4.  CD  framework  process  [2]. 


2.4.1  Pre-processing 

The  form  of  pre-processing  depends  on  the  data  under  analysis.  For  example,  if  the 
HS  image  data  is  “raw”  (obtained  directly  from  the  sensor),  then  atmospheric  com¬ 
pensation  and  band  reduction  (related  to  specific  atmospheric  distortion  windows) 
should  be  performed.  Typical  pre-processing  steps  include:  georegistration  ED  ,  or- 
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thorectification  [32],  atmospheric  compensation  {e.g.,  flat  fleld  correction,  empirical 
line  calibration,  QUick  Atmospheric  Correction  (QUAC)  and  Fast  Line-of-sight  At¬ 
mospheric  Analysis  of  Spectral  Hypercnbes  (FLAASH)  [33|).  Dimensionality  redac¬ 
tion  [331-137]  and  resolntion  enhancement  [38]. 

Georegistration  is  the  alignment  of  an  nnreferenced  image  with  geodetic  imagery. 
In  general,  this  process  aligns  two  images  and  correlates  the  images  to  a  physical 
location  by  examining  a  set  of  distingnishable  points.  If  not  performed  properly  it 
can  have  a  negative  effect  on  the  performance  of  the  CD  system.  For  example,  if 
the  accnracy  of  the  georegistration  degrades  snbst antially,  it  may  not  be  possible  to 
correlate  spatial  pixels  in  the  imagery  to  known  points  on  the  gronnd.  It  wonid  then 
become  nearly  impossible  to  obtain  an  accnrate  location  change  resnlt  referenced  to 
a  gronnd  position. 

Orthorectiflcation  is  the  process  of  registering  data  pixels  to  a  map  grid.  The 
end  resnlt  is  that  the  often  visnally-distorted  data  cnbes  captnred  from  a  HS  sensor 
are  spatially  re-sampled  so  that  they  may  be  overlaid  on  a  common  set  of  gronnd 
coordinates.  It  then  becomes  possible  to  measnre  distances  between  objects  in  the 
imagery  since  the  location  of  every  pixel  is  known.  The  drawback  to  orthorectiflcation 
is  that  the  spatial  interpolation  rontines  alter  the  spectral  content  of  the  HS  vectors 
within  the  image  cnbe  which  may  degrade  detection  performance.  Fig.  [^demonstrates 
a  notional  orthorectiflcation  scenario.  The  left  image  displays  an  object  as  seen  on  the 
gronnd,  known  to  be  a  sqnare.  The  middle  image  is  the  raw  image  with  the  object  as 
captnred  by  the  sensor  to  inclnde  geometric  distortion  (thns  the  sqnare  is  distorted). 
The  right  image  is  representative  after  orthorectiflcation;  the  image  is  transformed  so 
that  the  sensor  captnred  the  object  from  a  perpendicnlar  view. 
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Figure  5.  (left)  Square  object  as  it  appears  on  the  ground,  (middle)  square  object  sam¬ 
pled  with  a  pushbroom  sensor  (horizontal  distortion),  (right)  resulting  orthorectified 
imagery. 

Atmospheric  scattering,  absorption,  and  reflections  as  well  as  local  environmental 
reflections  contaminate  the  radiance  received  at  the  sensor  (Fig.  |^.  The  interpre¬ 
tation  of  surface  features  is  enhanced  by  correcting  for  these  effects.  Atmospheric 
compensation  is  a  technique  used  to  remove  the  effects  of  the  gases,  aerosols,  and 
water  vapor  present  in  the  atmosphere  on  the  spectra  observed  by  the  HS  sensor. 
This  is  particularly  important  for  interpreting  temporal  changes  in  imagery  taken  at 
a  different  time  of  day  or  a  different  time  of  year.  For  a  visible  to  near  infrared/short 
wave  infrared  (VNIR/SWIR)  HS  sensor,  there  are  five  intervening  atmosphere/envi¬ 
ronment  effects  to  the  spectral  radiance  that  are  ultimately  received  by  the  sensor  [39] : 

1.  Direct  solar  irradiance  reflected  off  the  object  and  transmitted  through  the 
atmosphere  to  the  sensor; 

2.  Indirect  solar  irradiance  scattered  by  the  atmosphere,  reflected  off  the  object, 
and  transmitted  through  the  atmosphere  to  the  sensor; 

3.  Scattered  irradiance  from  the  ground  and  local  objects  reflected  off  the  object 
and  transmitted  through  the  atmosphere  to  the  sensor; 

4.  Ground  scattered  radiance  scattered  further  by  the  atmosphere,  reflected  off 
the  object,  and  transmitted  through  the  atmosphere  to  the  sensor,  also  known 
as  adjacency  effect;  and 
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5.  Upwelling  path  radiance  scattered  by  the  atmosphere  directly  to  the  sensor. 


The  effect  of  the  atmosphere  on  the  solar  irradiance  spectrum  is  illustrated  in  Fig. 
Note  the  relative  loss  in  intensity  at  several  wavelengths  between  900  nm  and  1100  nm 
due  to  gas  and  vapor  absorption.  Also  note  that  two  major  absorption  features 
centered  around  1400  nm  and  1900  nm  are  predominantly  due  to  atmospheric  water 
vapor. 


Figure  6.  Solar  reflective  effects  described  by  (1)  Direct  solar  reflectance,  (2)  Indirect 
solar  reflectance,  (3)  Local  object  scattering,  (4)  Adjacency  effect,  and  (5)  Scattered 
path  radiance. 
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Figure  7.  Solar  irradiation  curves:  (dashed)  blackbody  radiator  at  5900K  (estimate  of 
the  suns  solar  radiation),  (solid  attenuated)  solar  irradiation  curve  at  sea  level,  (solid 
non-attenuated)  solar  irradiation  curve  at  the  top  of  the  atmosphere  (energy  from  the 
sun  entering  the  atmosphere)  [3]. 

Dimensionality  reduction  plays  a  role  in  reducing  processing  time  [Ml  EZ],  im¬ 
proving  visualization  [35l  [36] ,  and  improving  accuracy  results  in  the  assessment  |S7|. 
It  is  often  desirable  to  remove  some  of  the  bands  from  the  HS  data  cube  that  have 
low  Signal-to-Noise  Ratio  (SNR);  the  bands  corresponding  to  atmospheric  water  ab¬ 
sorption  are  commonly  among  those  with  low  SNRs  as  the  gases  and  vapor  in  the 
atmosphere  between  the  (airborne)  sensor  and  the  ground  plane  tend  to  absorb  light 
at  those  wavelengths. 

2.4.2  Change  detectors 
2.4. 2.1  Overview 

A  comprehensive  survey  of  image  CD  algorithms  is  presented  in  |6|.  A  variety 
of  change  detectors  have  been  developed,  and  many  have  been  summarized  and  case 
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studied  [T^l  ESI  EOHIS].  Many  of  the  documented  algorithms  can  be  modihed  for 
use  with  HS  images,  though  few  actually  develop  the  detector  specihcally  for  HSCD. 
Deer  |13]  categorizes  CD  methods  based  on  the  notion  of  pixel-,  feature-,  and  object- 
level  CD.  With  the  advent  of  HS  images  a  fourth  category  of  change  detectors  has 
emerged,  namely  spectral-based  change  detectors.  A  summary  of  these  four  categories 
of  CD  are  listed  in  Table  |2l 

Table  2.  Change  detectors. 


Pixel-based 

Feature-based 

Object-based 

Spectral-based 

—Image 
differencing 
—Image  rationing 
—Image  regression 
—Change  vector 
analysis  (CVA) 

— Endmember 
analysis 

—Principal 
components 
analysis  (PC A) 
—Local  texture 
—Shape  analysis 
—Vegetation  index 
differencing 
(NDVI)  -Wavelet 

—Artificial 
intelligence 
—Artificial  neural 

networks 

—Direct  multidate 

classification 

—Fuzzy 

post-classification 
comparison 
— Post-classification 
comparison 

—Spectral  angle 

mapper 

—Spectral 

correlation  mapper 

—Model-based 

chronochrome 

—Model-based 

covariance 

equalization 

Pixel-level  CD  requires  additional  processes  as  the  computations  involve  numeric 
values  of  each  image  band  {e.g.,  taking  image  differences  or  performing  image  ratios). 
These  techniques  evaluate  only  the  difference  between  corresponding  pixels  between 
the  reference  and  test  images.  Pixel-level  CD  requires  less  computation  since  only  one 
pixel  is  considered  at  a  time.  However,  it  is  very  sensitive  to  noise  and  illumination 
differences  since  it  does  not  exploit  local  structure  information. 

Feature-level  CD  involves  a  transformation  of  the  spectral  or  spatial  properties  of 
the  image  {e.g.,  principal  component  analysis,  texture  analysis,  or  vegetation  analy¬ 
sis). 

Object-level  CD  requires  the  most  advanced  processing  {e.g.,  post-classihcation 
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processing  or  neural  networks  analysis).  The  use  of  machine-learning  algorithms, 
artihcial  neural  networks,  and  decision  tree  classihers  have  gained  attention  in  recent 
years,  has  become  an  alternative  to  conventional  CD  approaches  [lU  Il5] .  Increased 
classihcation  accuracy  is  often  cited  as  the  primary  reason  for  developing  and  applying 
these  techniques;  however,  these  approaches  are  computationally  complex  and  can 
require  a  considerable  number  of  trusted  training  samples.  Advantages  of  object- 
level  CD  routines  include:  no  need  for  post-CD  hltering  and  smoothing;  less  impact 
of  slight  geometric  offsets  between  image  datasets;  and  the  ability  to  include  context 
relationships  to  improve  CD  results  |1U].  The  disadvantages  are  that  the  processes 
are  laborious  and  usually  require  previous  knowledge  of  the  study  area  as  well  as 
accurate  dehnition  of  the  changes  to  be  identihed. 

Spectral-based  detectors  make  a  direct  comparison  between  spectra  utilizing  the 
entire  spectrum  of  information  to  provide  areas  of  change.  There  exists  the  potential 
to  implement  a  spectral  matching  change  detector  using  techniques  typically  formu¬ 
lated  for  spectral  detection  m-  Here  one  exchanges  the  spectral  library  for  a  spectra 
signature  at  each  corresponding  spatial  pixel  in  the  test  image.  However,  these  meth¬ 
ods  inherently  assume  perfect  registration  since  they  perform  pixel-level  CD. 

Junior  et  al.  |18]  explore  the  spectral  angle  mapper  (SAM)  and  spectral  correlation 
mapper  (SCM)  classihcation  methods  and  use  them  as  spectral  change  detectors 
(SCD).  SCD,  such  as  the  Chronochrome  (CC)  and  Covariance  Equalization  (CE)  |19] 
algorithms,  perform  model-based  CD  through  signature  prediction. 

2.4. 2. 2  Chronochrome  change  detector 

The  chronochrome  change  detector  is  a  fundamental  component  of  this  disserta¬ 
tion  and  a  derivation  of  that  method  is  provided  here.  In  general,  the  chronochrome 
provides  the  linear  least  mean  squares  solution  to  the  prediction  error,  commonly 
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known  as  the  optimal  Wiener  filter  solntion  [SU].  The  chronochrome  reqnires  the 
cross-covariance  between  the  image  pair  because  it  is  a  pixel-to-pixel  computation 
forcing  precise  registration  which  is  difficult  to  achieve.  It  does,  however,  use  the 
generalized  likelihood  ratio  test  (GLRT),  which  is  well  understood  in  the  detection 
community,  and  unknown  parameters  are  replaced  with  their  maximum  likelihood  es¬ 
timates.  Covariance  equalization  is  an  alternative  to  chronochrome  by  approximating 
the  cross-covariance  through  whitening  in  order  to  eliminate  the  cross-statistics  [51]. 

Under  similar  environment  conditions  and  perfect  registration  where  no  parallax 
exists  then  the  CD  hypothesis  test  is 

Ho  :  t  =  r. 

Hi  :  t  ^  r, 

where  Hq  is  the  null  hypothesis  of  no  change.  Hi  is  the  alternative  hypothesis  of  a 
change,  and  r  and  t  are  two  measurements  from  a  reference  image  and  test  image. 
In  the  formulation,  lower-case  boldface  characters  represent  vectors,  and  upper-case 
boldface  characters  represent  matrices.  However,  in  reality,  the  respective  signatures 
can  be  slightly  different  as  a  result  of  intrinsic  variability.  As  a  result,  the  measure¬ 
ment  for  a  given  pixel  at  time-one  (r)  may  be  signihcantly  different  than  that  at 
time- two  (t),  even  when  the  underlying  material  associated  with  the  pixel  has  not 
changed.  Thereby,  the  measurements  are  assumed  to  have  a  linear  transformation, 
and  thus,  the  hypothesis  test  is 

Ho  :  t  =  Ar  -f  b  -|-  n. 

Hi  :  t  ^  Ar  -f  b  -|-  n, 

where  A  is  the  unknown  gain  matrix,  b  is  an  unknown  offset  vector,  and  n  is  additive 
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Gaussian  noise.  The  unknown  parameters  are  assumed  space-invariant  such  that 
the  global  sample  statistics  of  r  and  t  are  sufficient  to  estimate  the  transformation 
parameters. 

The  estimated  gain  matrix  A  and  offset  vector  b  are  determined  from  a  mini¬ 
mum  mean-squared  error  perspective  following  the  standard  Wiener  hlter  derivation 
requiring  the  error  to  be  orthogonal  to  the  reference  spectra  r,  or 

E  [[t  —  (Ar  -I-  b)]r'^}  =  0.  (1) 

Additionally,  the  residual  error  is  required  to  be  zero  mean, 

J5{[t-(Ar  +  b)|}=0,  (2) 

The  estimate  for  b  is  hrst  found  in  terms  of  A  as 

E  {[t  -  (Ar  b)]}  =  0 
E[t]  -  AE[t]  -  b  =  0 
b  =  E[t]  -  AE[r] 

b  =  lilt  —  Ariir  (3) 

where  liir  and  liit  are  the  sample  mean  vectors  from  time-one  and  time-two.  The 
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estimate  of  A  is 


E  |tr^  —  Arr^  —  br^}  = 

F£'[tr'^]  —  Ai?[rr'^]  —  bi?[r'^]  = 

£'[tr'^]  —  Aii^[rr'^]  —  (m^  —  Amr)-E'[r'^]  = 

£'[tr'^]  —  Aii^[rr'^]  —  mt-E'[r'^]  +  Amr-E[r'^]  = 

£^[tr'^]  —  mt-E[r'^]  = 

£^[tr'^]  —  mt-E[r'^]  = 

A  = 

where  Cr  is  the  sample  covariance  matrix  of  the  time-one  data  and  Ctr  is  the  sample 
cross-covariance  of  the  time-one  and  time-two  data. 

After  applying  the  transformation  to  the  test  image  nsing  the  estimates  derived, 
an  error  matrix  is  formed  via  differencing  the  reference  and  prediction  images  155. 

In  order  to  address  registration  concerns  and  create  an  algorithm  that  is  more 
robust  to  registration  errors,  the  covariance  equalization  algorithm  was  developed 
Covariance  equalization,  which  is  based  on  the  whitening  principle,  was  introduced 
as  an  approximation  to  the  chronochrome  algorithm,  and  it  has  a  number  of  practical 
advantages  in  terms  of  implementation. 

To  derive  the  covariance  equalization  algorithm,  follow  chronochrome  derivation 
and  then  begin  with  Eq.  Q,  restated  here. 


The  cross-covariance  can  be  approximated  as  the  covariance  of  Ct  [52].  Each  image 
is  brought  into  a  new  whitened  space  based  on  their  own  covariance  matrix.  By  doing 


0 

0 


AE'[rr'^]  —  AmrE[r'^] 

A(E[rr'^]  —  mrE'[r'^]) 
E[tr'^]  -  mtm^  _  ^ 
Elrr'^l  —  mriii^  Cr 


(4) 
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this,  registration  is  less  of  a  concern  in  the  prediction  stage  because  it  does  not  require 
the  cross-covariance. 

Eigenvalue  decomposition  refers  to  the  factoring  of  a  matrix  into  a  diagonal  ma¬ 
trix  of  eigenvalues,  and  corresponding  eigenvectors.  Since  Ct  and  Cr  are  covariance 
matrices  and  by  dehnition  are  square,  symmetric,  and  positive  semi-dehnite,  they  can 
be  written  as. 


Ct  =  = 

Cr  =  $rArd>^  =  , 

where  and  <hr  are  orthogonal  eigenvector  matrices  and  A*  and  A^  are  diagonal 
matrices  of  eigenvalues  resulting  from  the  eigenvalue  decomposition  of  the  covariance 
matrices.  Then,  the  symmetric  square  root  of  the  covariance  matrices  can  be  whitened 
through  the  transformation  <hA“^A<|)T’ 


Ct  = 

Cr  = 

$rAy"C- 

The  estimate  for  covariance  equalization  of  A  is  then 

A  = 

Ct 

Cr 

(5) 

(6) 

A  = 

cX 

A  = 

cX 

(7) 

A  = 

cXcX'^ 

(8) 

To  equalize  the  mean  vectors,  the  estimate  for  b  is  given  again  by  Eq.  (|^. 
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2. 4. 2. 3  Local  co-registration  adjustment 


To  account  for  co-registration  errors  between  the  reference  and  test  images,  Theiler 
and  Wohlberg  [53]  describe  a  method  that  extends  any  test  statistic  for  anomalous 
change  and  minimizes  that  test  statistic  for  pixels  in  the  reference  image  that  are  in  the 
neighborhood  of  the  test  image  pixel  under  evaluation.  It  does  so  by  comparing  the 
measurement  of  vector  r  in  the  reference  image  with  a  3x3  window  of  its  correspond¬ 
ing  vector  t  in  the  test  image.  Then,  the  nearest  neighbor  with  the  most  anomalous 
relationship  is  the  mis-registered  compensated  pixel.  Local  Co-Registration  Adjust¬ 
ment  (LCRA)  includes  both  symmetric  (independent  of  the  direction  of  minimization) 
and  sub-pixel  variants  of  the  algorithm.  It  assumes  that  a  true  anomalous  change 
is  independent  of  the  direction  of  the  minimization  and  therefore  the  minimization 
is  performed  in  both  directions  {i.e.,  r  to  t,  or  vice  versa)  where  the  maximum  of 
the  resulting  minima  is  selected  as  the  hnal  results.  Finally,  the  sub-pixel  LCRA 
redehnes  the  optimization  windows  to  include  fractional  pixel  shifts.  This  is  applied 
using  interpolation  to  resample  the  images  at  higher  resolution  prior  to  processing. 

2.4. 2.4  Change  detector  summary 

The  notions  of  changes  that  are  signihcantly  different  or  unimportant  vary  by 
application,  which  makes  it  difficult  to  directly  select  one  from  the  plurality  of  CD 
routines.  There  is  no  optimal  CD  routine  and  the  procedure  that  is  most  appropriate 
depends  upon  the  application  (persistent  ISR,  target  tracking,  etc.)  and  the  data 
(image  quality,  geography  of  study  area,  etc.).  In  addition,  the  choice  of  CD  routine 
for  the  application  should  include  a  priori  knowledge  of  the  OCs  in  order  to  dehne 
the  appropriate  techniques  that  reveal  signihcant  change.  This  aspect  is  important 
for  CD  as  it  eventually  determines  if  the  data,  the  features,  or  the  detector  are  cost 
effective. 
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2.4.3  Thresholding 


Thresholding  is  a  fundamental  method  applied  in  many  image  processing  and 
detection  methods.  If  the  change  output  level  exceeds  a  given  threshold,  the  pixel  is 
regarded  as  a  change  pixel.  Setting  the  threshold  is  a  crucial  step.  If  set  too  high, 
it  will  suppress  the  detection  of  signihcant  changes.  If  set  too  low  it  will  swamp 
the  change  map  with  spurious  changes.  A  typical  thresholding  routine  is  modeled  in 
Fig.|8} 


Figure  8.  Typical  thresholding  routine. 


Per  the  framework  diagram  in  Fig.  responses  from  the  detector  optionally  un¬ 
dergo  a  thresholding  stage  directly  or  after  post  processing.  When  accomplished 
directly  after  the  detector,  it  serves  to  reduce  the  amount  of  information  processed  in 
the  post-processing  stage.  If  accomplished  after  the  hrst  stage  of  post-processing,  it  is 
used  to  rehne  a  hltered  response  (due,  e.g.,  to  classihcation)  for  potentially  improved 
results.  Regardless  of  when  the  threshold  is  applied,  before  or  after  post-processing, 
the  method  produces  a  binary  change  mask.  As  such,  potential  objects  of  interest 
(signihcant  changes)  are  separated  from  the  background.  The  binary  image  mask 
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{B)  identifies  the  changed  region  [5j  where 


{1  if  signihcant  change, 

0  otherwise. 

For  most  CD  methods  in  the  literature,  the  threshold  routine  is  not  explicitly 
dehned  and  there  is  no  theoretical  way  for  setting  it  ahead  of  time  because  of  the 
complexity  of  the  background  clutter.  Thresholds  are  typically  computed  globally 
or  locally  and  often  based  on  hypothesis  testing.  Global  threshold  routines  operate 
on  all  pixels  in  the  image  whereas  local  threshold  routines  operate  on  a  subset  of 
the  change  information.  Both  approaches  use  spatial  and/or  spectral  information  for 
the  assessment.  However,  the  spectral  information  tends  to  be  more  robust  in  high 
noise  environments.  Ortiz-Rivera  et  al.  la  demonstrated  that  the  local  thresholding 
methodology  improved  CD  results  in  their  HSCD  experiment. 

Knowing  what  to  look  for  in  an  image  might  simplify  estimating  the  threshold, 
but  trial-and-error  is  the  most  common  approach.  Several  thresholding  methodolo¬ 
gies  have  been  proposed  in  the  literature;  however,  few  of  them  are  specihc  to  CD. 
Rosin  [21]  describes  four  different  methods  for  selecting  CD  thresholds;  either  the 
noise  or  the  signal  is  modeled,  and  the  model  covers  either  the  spatial  or  intensity 
distribution  characteristics.  Spatial  features,  such  as  texture  and  local  variance,  can 
help  with  establishing  the  CD  threshold  [5l|,  but  this  implies  that  user  objectives  are 
determined  prior  to  the  design  of  the  change  analysis  problem.  Rosin  does  conclude 
that  the  most  promising  methods  are  spatial  but  adds  more  extensive  testing  and 
quantitative  assessment  is  required. 

Regardless  of  the  specihc  technique  used,  a  priori  knowledge  regarding  the  data 
and  the  application  is  helpful  in  selecting  a  relevant  threshold.  Approaches  to  de¬ 
termine  the  threshold  values  are  still  an  active  area  of  research  and  are  most  often 
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integrated  into  the  post-processing  stage. 


2.4.4  Post-processing 

Post-processing  can  occur  on  the  change  output  prior  to  thresholding  or  after 
thresholding  on  the  hltered  change  responses  (Fig.  [^.  Post-processing  accomplished 
at  either  point  serves  to  improve  the  candidate  change  mask.  The  initial  binary  mask 
generated  after  thresholding  can  highlight  insignihcant  changes  and  post-processing 
rehnes  the  results  potentially  improving  accuracy. 

The  simplest  techniques  post-process  the  change  mask  with  standard  binary  image 
processing  operations.  Morphological  operations  are  commonly  used  as  a  tool  for 
extracting  information  from  an  image  to  describe  regional  shapes;  morphology  is  a 
way  of  giving  meaning  to  a  certain  region  in  an  image  [B].  Post-processing  applied 
to  the  binary  image  has  the  advantage  of  reducing  false  alarm  probabilities  at  low 
computational  cost.  However,  standard  binary  image  processing  operations  may  not 
always  be  valid  in  CD  |55j. 

Zhang  et  al.  |23]  use  the  minimum  noise  fraction  (MNF)  to  post-processing  the 
raw  CD  response  in  order  to  separate  signal  from  noise.  The  change  output  is  then 
processed  with  a  Markov  Random  Field  (MRF)  model  [56]  for  the  threshold  step  in¬ 
corporating  contextual  information  to  improve  accuracy  and  reliability.  Lee  et  al.  123 
investigate  the  feasibility  of  incorporating  symbolic  reasoning  to  remove  insignihcant 
changes  due  to  shadows,  clouds  and  partial  occlusion  of  existing  objects.  Their  over¬ 
all  CD  system  concept  includes  a  change  interpretation  node  which  applies  rules  in 
a  hypothesis-driven  fashion  to  the  change  output  to  determine  the  relevance  of  that 
change.  In  general,  the  symbolic  reasoning  concept  is  not  well  dehned,  nor  does  it 
represent  a  variety  of  applications. 

Decision  fusion  is  also  used  in  post-processing  to  mitigate  false  alarms.  Decision 
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fusion  is  the  combination  of  two  or  more  algorithms  and  is  most  often  applied  to  the 
change  detector.  The  results  of  each  CD  routine  are  compared  and  their  contribution 
to  the  hnal  decision  is  weighted  to  be  more  or  less  signihcant  than  the  others.  A 
decision  fusion  step  logically  combines  the  hnal  change/no  change  locations  then  cor¬ 
relates  those  from  each  change  detector.  All  weighted  outputs  are  then  thresholded, 
at  which  point  additional  false  alarms,  and  perhaps  some  true  detections,  are  hltered 
out.  Fig.  provides  an  illustration  of  a  conventional  decision  fusion  system.  For 
each  pixel  in  the  change  image  there  is  consensus  as  to  whether  or  not  it  contains 
signihcant  change  from  the  reference  image.  (This  concept  and  hgure  was  motivated 
by  Hall  and  Llinas’s  |1]  tutorial  on  data  fusion.) 


Figure  9.  Conventional  decision  fusion  schema.  Image  motivated  by  Hall  and  Lli¬ 
nas’s  [4]. 

Other  post-processing  false  alarm  mitigation  routines  use  geographic  information 
system  (GIS)  data  [571  EH]  to  enable  the  delivery  of  quantihed  or  stratihed  change 
maps  that  are  consistent  with  change  product  delivery  |in|.  (A  GIS  is  a  system  that 
captures,  stores,  analyzes,  manages,  and  presents  data  that  are  linked  to  location.) 
The  advantage  of  using  GIS  is  the  ability  to  incorporate  diherent  data  sources  into 
change  analysis  methods,  thereby  providing  a  better  understanding  and  identihcation 
of  the  changes  that  have  occurred  in  a  geographic  area  of  interest.  The  more  promi¬ 
nent  work  by  Zhang  et  al.  [5H|  puts  forward  a  holistic  strategy  of  GD  based  on  remote 


sensed  images  and  GIS  data  and  claims  it  addresses  the  shortcomings  of  traditional 
methods,  such  as  the  error  accumulation  and  error  transfer  caused  by  pre-processing 
techniques  such  as  data  co-registration  and  feature  extraction.  The  result  is  improved 
accuracy  and  reliability  of  the  change  results. 

Post-processing  can  also  be  in  the  form  of  visual- quality  improvement  to  aid  the 
image  analyst  in  assessing  CD  results.  Much  of  the  hltering  and  post-processing  is  per¬ 
formed  to  enhance  visual  characteristics  of  images.  HS  visualization  is  a  distinct  but 
related  problem  to  band  reduction.  Visualization  techniques  render  either  enhanced 
true  color  or  false  color  images  derived  from  the  HS  image,  typically  in  an  unsuper¬ 
vised  manner  [36].  Resolution  enhancement  aims  to  compensate  for  the  typical  low 
spatial  resolution  of  HS  sensors,  primarily  by  fusing  additionally  available  imagery 
or  performing  spatial  superresolution  [59].  The  result  is  a  HS  cube  with  increased 
spatial  resolution  which  greatly  aids  visualization  for  inquiry  of  the  performance  in 
CD  results. 

2.4.5  Analysis/ Assessment 
2.4. 5.1  Overview 

Analysis  and  assessment  is  important  when  developing  a  CD  system  and  that 
assessment  occurs  by  comparing  results  with  a  truth  image  mask  (a  verihed  ground 
truth  pixel  map).  Accuracy  assessment  techniques  in  CD  originate  from  those  of 
remote  sensing  image  classihcation.  It  is  natural  to  extend  the  accuracy  assessment 
techniques  for  processing  a  temporarily  stationary  image  to  that  of  multi-temporal 
images. 

Rosin  [60]  comments  that  accurate  ground  truth  is  essential  either  through  the 
use  of  synthetic  imagery  or  annotation  and  an  image  analyst.  However,  it  is  often 
difficult,  if  not  impossible,  to  collect  reliable  ground  reference  due  to  cost,  time. 
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and  accessibility  of  the  study  site.  Three  general  approaches  to  obtain  ground  truth 
references  are  held  survey,  acquisition  of  airborne  high-resolution  imagery  (HRI), 
and  visual  interpretation.  The  hrst  two  approaches  incur  a  large  cost  and  the  third 
may  not  be  accurate.  Since  the  ground  truth  data  provides  a  reference  point,  the 
thresholding  aspect  in  CD  can  be  tested  for  correctness  of  the  change  results.  As 
described  previously,  incorporating  CIS  information  can  improve  the  end-to-end  CD 
system  results.  CIS  further  provides  convenient  tools  for  multi-source  data  processing 
and  is  effective  in  handling  the  CD  analysis  using  multi-source  data  0. 


2.4. 5. 2  Error  matrix,  receiver  operating  characteristic  (ROC)  curve, 
and  area  under  the  ROC  curve  (AUC) 

The  performance  of  a  CD  algorithm  can  be  evaluated  visually  and  quantitatively 
based  on  application  needs.  The  truth  reference  and  a  visual  interpretation  com¬ 
bination  can  provide  a  sense  of  performance.  However,  an  image  analyst  may  hnd 
that  the  change  map  requires  extensive  examination  of  imagery  to  assess  the  relation 
of  the  true  change.  Among  the  quantitative  assessment  techniques  for  analysis  of 
change  masks,  the  most  efficient  and  widely  used  are  the  error  matrix  [13],  the  re¬ 
ceiver  operating  characteristic  (ROC)  [61]  curve,  and  the  area  under  the  ROC  curve 
(AUC)  [62]. 

Oort  [33]  describes  the  common  change/no  change  error  matrix  used  in  analyzing 
results.  Here,  CD  is  equivalent  to  a  two-class  classihcation  problem.  A  typical  error 


matrix  for  a  two-class  classihcation  problem  is  illustrated  in  Fig.  ]R  As  shown,  the 
error  matrix  is  compiled  into  a  table  with  two  rows  and  two  columns  that  report  the 
number  of  false  positives,  false  negatives,  as  well  as  the  complementary  true  positives 
and  true  negatives. 
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Classes  assigned  by  change 
detector/threshold 


True 

Positives 

(TP) 

False 

Positives 

(FP) 

TP+  FP 

False 

Negatives 

(FN) 

True 

Negatives 

(TN) 

FN  +TN 

TP+FN 

FP  +  TN 

TP  +  TN 

Figure  10.  Typical  change  detection  error  matrix. 


The  error  matrix  provides  a  means  to  calculate  accuracy  from  both  the  perspective 
of  the  system  and  the  consumer.  The  accuracy  rate  refers  to  the  percentage  of  correct 
predictions  when  compared  with  the  actual  classifications  in  the  test  data.  Addition¬ 
ally  accuracies,  such  as  producer’s  and  consumer’s  accuracy,  are  easily  computed  from 
the  error  matrix.  Producer’s  accuracy  is  an  expression  of  errors  of  omission,  or  false 
dismissals  (for  example,  from  the  perspective  of  CD,  labeling  a  pixel  classified  as  no 
change  when  in  fact  it  is  change).  Consumer’s  accuracy  is  an  expression  of  errors  of 
commission,  false  alarms  (for  example,  again  from  the  perspective  of  CD,  labeling  a 
pixel  as  change  when  in  fact  it  is  not). 


ROC’s  (Fig.  11)  are  generally  interpreted  in  one  of  two  ways.  If  ROC  one  domi¬ 
nates  ROC  two,  then  system  one  performs  better  than  system  two  across  the  entire 
operating  range  [6H  |65].  If  ROC  one  and  ROC  two  cross,  the  interpretation  is  that 
system  one  performs  better  than  system  two  over  different  portions  of  the  operating 
range  |66].  ROC  curves  are  obtained  from  binary  hypothesis  testing  by  using  several 
threshold  values  [71E7]. 

Bradley  in  [02]  (and  others)  suggest  using  the  area  under  the  ROC  curve  (AUC). 
The  ROC  is  generated  in  the  normal  way  and  the  area  is  estimated  using  techniques 
such  as  the  trapezoidal  integration  method  |62|.  The  AUC  is  not  without  its  problems. 
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Namely,  it  summarizes  the  ROC  curve  with  a  single  number  and  loses  the  ROC’s 
granularity.  It  does,  however,  provide  a  way  to  distinguish  between  two  systems 
when  the  ROC  of  one  system  crosses  the  ROC  of  another  system  (he.,  one  ROC  does 
not  dominate  the  other  ROC). 


ROC  Curves 


Probability  of  False  Alarm 

Figure  11.  Three  hypothetical  ROC  curves.  The  probability  of  detection  (Pd)  is  plotted 
against  the  probability  of  false  alarm  (Pfa)  based  on  changes  to  a  detection  threshold. 
A  value  of  0.5  (blue  line)  suggests  random  guesses  while  a  value  of  1  indicates  correct 
classification;  as  the  curves  (red  and  green)  approach  the  value  of  Pp  =  1  and  Pfa  —  0, 
the  detector  performance  is  said  to  improve. 


In  order  to  evaluate  the  quality  of  a  change  mask  independent  of  the  choice  of 
thresholding,  the  evolution  of  the  detection  probability  as  a  function  of  the  false 
alarm  probability  may  be  evaluated  against  the  ground  truth  data.  The  ground 
truth  is  scored  against  the  binary  threshold  image  and  can  include  different  scoring 
methodologies  such  as  proximity,  touch,  target  centroid,  blob  centroid,  vote  and  pixel- 
level  (Fig.  [I^. 
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Proximity 

♦  A  blob  within  the  immediate  vicinity  of  the  target  counts  as  a  detection. 

♦  This  method  only  requires  the  center  plus  radius  target  mask. 

♦  It  is  good  for  drawing  attention  to  a  region  containing  targets. 

♦  This  allows  large,  sprawling  blobs  to  be  counted  as  detections. 


Touch 

♦  The  blob  of  pixels  must  touch  the  target. 

♦  This  is  the  Proximity  Method  with  tight  boundaries. 

♦  The  target  masks  must  be  regular  or  irregular  polygons. 

♦  This  method  allows  large,  sprawling  blobs  to  be  counted  as  detections. 


Target  Centroid 

♦  The  blob  must  touch  the  center  of  target. 

♦  It  ensures  that  some  pixels  are  on  the  target. 

♦  This  method  only  needs  the  center  position  of  the  target  mask. 

♦  This  allows  large,  sprawling  blobs  to  be  counted  as  detections. 


Blob  Centroid 

•  The  center  of  the  blob  must  touch  the  target. 

♦  This  is  one  of  the  most  conservative  approaches. 

•  It  ensures  that  some  pixels  are  on  the  targets. 

♦  This  method  works  well  with  any  type  of  target  mask. 


Vote 

♦  There  must  be  a  certain  percent  of  pixels  on  the  target  to  count  as  a 
detection. 

♦  This  method  should  have  detailed  target  masks  (irregular  polygons). 

♦  This  procedure  must  evaluate  performance  on  a  pixel-by-pixel  basis. 


Pixel  Level 

♦  Any  pixel  on  the  target  mask  counts  as  a  detection. 

♦  This  is  the  Vote  Method  without  a  percent-of-pixels  threshold. 

♦  It  is  extremely  lenient;  only  one  pixel  is  required  to  touch  the  target. 


Figure  12.  Possible  scoring  methodologies  used  during  assessment  [5]. 


Accuracy  assessment  provides  benefits  in  the  CD  framework  including  the  decision 
of  choosing  the  specihc  methodologies  of  each  segment  of  the  processing  chain  and 
reporting  the  overall  reliability  of  the  CD  system.  There  is  still  a  need  for  research  in 
accuracy  assessment  for  CD,  especially  where  the  ground  reference  data  is  insufficient 
or  impossible  to  obtain. 

2.4.6  Discussion 

The  contribution  of  the  framework  presented  here  is  another  step  toward  under¬ 
standing  the  design  space  of  CD  systems  by  extracting  the  crucial  processing  steps 
in  each  processing  stage.  From  both  a  research  and  implementation  perspective,  the 
analysis  shows  how  each  stage  can  be  decomposed  and  implemented  as  part  of  an 
end-to-end  CD  system.  For  the  researcher,  the  framework  enables  one  to  understand 
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and  compare  their  individual  CD  research  efforts  (within  a  specific  method  specified 
in  the  framework)  in  different  end-to-end  CD  systems.  For  the  user,  the  framework 
helps  the  analyst  choose  specific  methods  for  each  process  in  the  CD  processing  chain 
to  achieve  best  results. 

2.5  Stereo  geometry 

Stereo  vision  is  often  considered  a  form  of  image  registration  because  the  aim  is  to 
correlate  a  secondary  image  (change  image)  with  a  reference  image  and  both  require 
some  form  of  correspondence  algorithm.  The  distinction  is  that  in  the  image  regis¬ 
tration  case,  the  process  lies  in  estimating  a  transformation  that  spatially  aligns  two 
images.  However,  in  stereo  vision  the  aim  is  to  utilize  the  stereo  geometry  to  locate 
the  disparity  between  corresponding  points,  allowing  the  creation  of  3-Dimensional 
(3D)  data.  Similar  to  image  registration,  stereo  vision  encounters  its  fair  share  of  diffi¬ 
culties  during  correspondence  matching.  Some  problems  commonly  endured  by  stereo 
algorithms  include  occlusions,  perspective  distortion,  and  illumination  variation. 

Stereo  vision  refers  to  the  ability  to  infer  information  on  the  3D  structure  and 
distance  of  a  scene  from  two  or  more  images  taken  from  different  viewpoints.  Typical 
stereo  vision  systems  solve  two  problems:  correspondence  and  reconstruction  [68]. 
Stereo  correspondence  is  a  well-understood  problem  in  computer  vision,  with  numer¬ 
ous  applications  including  the  previous  mentioned  image  registration  problem.  The 
reconstruction  problem,  given  a  number  of  corresponding  points  in  both  images  and 
geometry  of  the  stereo  pair,  permits  development  of  a  3D  map  of  the  scene. 

Stereo  matching  methods  try  to  solve  the  stereo  correspondence  problem.  This 
dissertation  aims  to  use  stereo  matching  to  identify  all  homologous  pixels  in  the 
reference  and  change  image  pair  to  yield  a  parallax  approximation  [691  ED].  The 
stereo  correspondence  problem  is  solved  by  matching  points  between  a  stereo  pair 
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of  images.  For  every  point  in  the  reference  image  of  the  stereo  pair,  a  matching 
point  is  found  in  the  test  image  of  the  stereo  pair.  This  match  is  described  by  a 
translation  vector  from  the  reference  pixel  in  the  reference  image  to  the  test  pixel  in 


the  test  image  (Fig.  13).  Typical  depth  applications  devise  a  disparity  map  from  the 
translation  vectors. 


Figure  13.  Depiction  of  residual  displacement  or  parallax. 


Matching  corresponding  points  in  each  image  can  be  a  computational  burden 
within  a  2-Dimensional  (2D)  array.  In  order  to  simplify  the  search,  epipolar  geometry 
is  often  exploited.  Epipolar  geometry  interprets  the  projective  relationship  for  the 
stereo  pair.  Hartley  nn  describes  the  epipolar  geometry  of  two  views  where  a  point 
in  one  view  dehnes  an  epipolar  line  in  the  other  view  on  which  the  corresponding 


point  lies.  Fig.  depicts  the  epipolar  geometry  for  a  stereo  pair.  The  3D  point 
p  is  imaged  in  two  views,  at  in  the  reference  image  and  p^  in  the  test  image. 
The  epipoles  are  designated  by  and  which  is  dehned  by  the  intersection  of  the 
line  joining  the  image  centers,  and  c*.  Thus  the  matching  of  points  is  somewhat 
simplihed  as  the  epipolar  geometry  enables  the  search  for  corresponding  points  only 
along  corresponding  image  lines.  The  epipolar  constraint  restricts  the  search  for  the 
match  of  a  point  in  one  image  along  the  corresponding  epipolar  line  in  another  image. 
Thus  the  search  for  correspondences  is  reduced  to  a  1-dimensional  (ID)  problem. 
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Alternatively,  false  matches  due  to  occlusions  can  be  rejected  by  verifying  whether  or 
not  a  candidate  match  lies  on  the  corresponding  epipolar  line. 


P 


imaged  in  two  views,  at  p,.  in  the  reference  image  and  pt  in  the  test  image.  The  epipoles 
are  designated  by  e^  and  e^  which  is  defined  by  the  intersection  of  the  line  joining  the 
image  centers,  Cr  and  Cj. 

Correspondence  algorithms  aim  to  describe  a  pixel,  or  area,  from  the  image  data 
that  can  be  used  to  efficiently  perform  a  similarity  measure  between  the  multi-view 
imagery.  There  are  two  steps  in  establishing  a  stereo  correspondence:  (1)  select 
(manually  or  automatically)  image  pixel  locations  to  match  and  (2)  determine  the 
similarity  of  that  match.  There  are  a  large  number  of  proposed  descriptors  for  auto¬ 
matic  pixel  selection  and  associated  similarity  measures  which  utilize  different  local 
image  features  {e.g.,  pixel  intensities,  color,  texture,  edges,  etc). 

In  performance  comparison  of  different  descriptors,  shows  that  (SIFT)  is  one 
of  the  most  stable  and  reliable  feature  detection  algorithms.  An  overview  of  the  SIFT 
keypoint  descriptor’s  construction  is  provided  where  a  complete  description  of  the 
SIFT  algorithm  can  be  found  in  [73].  The  SIFT  descriptor  is  based  on  a  monochro¬ 
matic  or  gray-scale  image  representation.  SIFT  features  are  local  histograms  of  gra¬ 
dient  directions  computed  over  different  parts  of  the  keypoint  region.  It  computes 
the  gradient  vector  for  each  pixel  in  the  keypoint ’s  pixel  neighborhood  and  builds 
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a  normalized  orientation  histogram  of  gradient  directions  to  create  a  128-dimension 
vector  for  feature  correspondence.  Orientation  invariance  is  achieved  by  estimating 
the  dominant  orientation  of  the  local  image  patch  and  is  normalized  for  invariance 
to  changes  in  illumination.  The  SIFT  process  of  generating  keypoints  is  performed 
in  two  steps:  (1)  feature  detection  and  (2)  feature  extraction.  Feature  detection 
or  keypoint  detection  is  simply  the  process  of  locating  salient  points  and/or  regions 
from  images,  in  order  to  construct  useful  local  image  descriptors.  Other  important 
attributes  of  a  reliable  and  meaningful  salient  keypoint  is  that  it  must  be  invariant 
to  image  transformation,  such  as  scale,  rotation,  and  affine  transform,  in  addition 
to  perspective  transformation,  illumination,  and  brightness  variations.  Feature  ex¬ 
traction  is  that  local  descriptors  are  computed  on  the  region  around  each  keypoint 
detected  by  the  keypoint  detection  process.  The  local  invariant  feature  descriptors 
are  formed  by  the  local  image  gradients  measured  at  the  selected  scale  in  the  re¬ 
gion  around  each  keypoint.  In  general,  extracted  features  must  be  highly  distinctive 
and  tolerant  to  image  noise,  changes  in  illumination,  uniform  scaling,  rotation,  and 
changes  in  viewing  direction. 

The  second  step  is  determining  the  match  through  a  similarity  measure  to  pro¬ 
vide  spatial  correspondence.  Various  similarity  measures  exist  {e.g.,  sum  of  absolute 
differences,  cross-correlation  coefficient,  mutual  information)  for  typical  feature  cor¬ 
respondence  methodologies.  The  objective  of  the  present  work  is  to  improve  HSCD 
by  eliminating  the  false  parallax  change.  As  such,  stereo  correspondence  is  utilized, 
not  for  image  registration  in  the  pre-processing  stage  of  the  CD  processing  chain,  but 
to  account  for  parallax  thus  eliminating  false  change.  Specihcally,  a  correspondence 
search  is  made  between  the  reference  and  test  image  pair  to  measure  the  parallax 
or  the  apparent  shift  of  a  3D  point,  p,  above  the  ground  plane  against  the  ground 
plane.  The  difference  between  the  position  of  the  pixel  in  the  reference  image  and  the 
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position  of  the  pixel  in  the  test  image  on  the  ground  plane  identihed  as  the  coordi¬ 
nates  of  the  corresponding  point  of  p  provide  an  estimate  of  the  parallax  shift.  The 
parallax  direction  associated  with  this  estimate  is  a  key  component  in  the  parallax- 
compensation  algorithm  developed  in  this  research. 
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III.  Change  Detection  Accounting  for  Mis- registration  and 

Parallax  Errors 


This  chapter  presents  a  new  method  to  account  for  image  mis-registration  and 
parallax-induced  false  detections  in  HSCD.  A  generalized  likelihood  ratio  test  (GLRT) 
is  constructed  to  provide  a  change  test  statistic  for  change  at  each  pixel  in  the  tempo¬ 
ral  hyperspectral  image  pair,  and  the  unit  square  restricted  form  of  bilinear  interpo¬ 
lation  d  is  used  to  account  for  estimates  of  image  mis-registration  for  the  sub-pixel 
case.  This  baseline  approach  is  then  modihed  to  incorporate  an  estimate  of  the  par¬ 
allax  direction  in  order  to  mitigate  false  alarms  due  to  parallax  effects  between  the 
same  two  hyperspectral  image  pairs.  The  organization  of  this  chapter  follows  this 
same  basic  construction. 

3.1  Accounting  for  image  mis-registration 

3.1.1  Deriving  the  mis-registration  generalized  likelihood  ratio  test 
(GLRT) 

In  the  following  derivations,  lower-case  boldface  characters  represent  vectors,  and 
upper-case  boldface  characters  represent  matrices.  Consider  a  hyperspectral  sensor 
that  makes  two  observations  of  a  scene  in  the  form  of  a  hyperspectral  image  composed 
of  A-element  spectra  made  at  spatial  locations  defined  by  an  image  coordinate  {x,y). 
The  image  coordinate  can  be  characterized  in  any  domain,  such  as  ground  plane 
or  angular  coordinates,  but  {x,  y)  is  used  to  represent  spatial  pixel  indices  in  this 
dissertation.  The  reference  image  is  denoted  as  v{x^y)  and  the  test  image  as  t{x,y) 
where  v{x,y)  G  3?^,  t{x,y)  G  3?'^.  Typically,  hyperspectral  images  are  indexed  as 
matrices  where  the  two-dimensional  spatial  nature  is  collapsed  into  one  dimension; 
however,  the  two-dimensional  spatial  relationships  are  explicitly  maintained  here  in 
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order  to  capture  spatial  mis-registration  impacts. 

Define  the  set  of  local  sub-pixel  mixtures  using  the  unit  square  restricted  form  of 
the  bilinear  mixture  model  around  a  reference  spectrum  0  v{x,y)  as 


{m[(a:,  y)]  Uy]  :  \aa:\  <  1,  \ay\  <  1}, 


(9) 


where 


m[(x,y);a^,aj  =  < 


r{x,y) 

(1  -  a^){l  -  ay)r{x,  y)  -h  (1  -  a:,)ayr{x,  y  +  1) 

-haa;(l  -  ay)r{x  -h  1,  y)  -f  a:,ayr{x  -h  1,  y  -h  f) 
(1  a^)(l  -  ay)Y{x,y)  (1  -h  a^)ayT(x,y  +  1) 

-a^(l  -  ay)T{x  -  1,  y)  -  a^ayT{x  -l,y  +  l) 
(1  -h  a^){l  +  ay)r{x,  y)  -  {1  +  a^)ayr{x,  y  -  1) 
-aa;(l  -h  a;y)r(x  -  1,  y)  -h  a^ayr{x  -  1,  y  -  1) 
(1  -  Q!^)(l  -h  Q!y)r(x,  y)  -  (1  -  a^)ayr{x,  y  -  1) 
-hax(l  +  ay)r{x  -h  1,  y)  -  a^ayr{x  +  l,y  -1) 


ii  ax  —  oiy  =  Q 


if  ax  >  0,ay  >  0 


if  Oa;  <  0,  Oy  >  0 


if  Oa;  <  0,  Oy  <  0 


if  Oa:  >  0,  Oy  <  0 


Assume  (ofa:,  ay)  are  unknown  but  can  range  from  zero  to  one  pixel.  Also  assume  that 
a  specihc  test  spectrum  t(a:,y)  maps  to  a  specihc  reference  spectrum  r(x,y)  within 
some  unknown  mis-registration  (Ax,  Ay). 

Consider  the  situation  where  there  is  some  unknown  mis-registration  (Ax,  Ay) 
between  the  reference  and  test  images  such  that  the  prediction  image  is  required  to 
be  of  the  form 

p(x,  y|  Ax,  Ay)  =  Ar(x  -f-  Ax,  y  +  Ay)  -|-  b.  (10) 


The  assumption  is  that  the  amount  of  mis-registration  can  vary  from  pixel  to  pixel 
according  to  some  known  or  assumed  statistics  dehned  by  a  probability  density  func¬ 
tion  (PDF)  /Ax,Ay(Ax,  Ay).  Letting  n  represent  the  sensor  noise,  the  hypothesis 
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test  corresponding  to  CD  is  written  as 


Ho  :  e{x,y)  =t{x,y) -p{x,y\Ax,Ay) +  n  =  0, 

Hi  :  e{x,y)  =t{x,y) -p{x,y\Ax,Ay) +  (11) 


The  GLRT  test  statistic  for  this  detection  problem  is  given  by  the  maximum  likelihood 
ratio  eg  dehned  as 


l{x,y) 


max  fe[e{x,y]A,h,Ax,Ay)\Hi] 
A,b,Ax,Aj/ 

max  fe[e{x,y]A,h,Ax,Ay)\Ho] 

A, b, Ax, At/ 


(12) 


where  /e  ( ■ )  is  the  probability  density  function  for  the  residual  error  (e)  in  Eq.  (12). 
Under  the  alternative  hypothesis,  for  anomalous  CD,  there  is  no  a  priori  information; 
therefore,  the  maximum  likelihood  estimates  are  the  measurements  themsleves  and 


the  numerator  in  Eq.  (12)  is  constant.  Eq.  (12)  is  written  as 


/e[e(a;,  r/;  A,  b.  Ax,  Ay)\Ho]  =  /e[e(x,  y;  A,  b)|ifo,  Ax,  Ay] /ax, av (Ax,  Ay).  (13) 


The  assumption  is  n  ~  A(0,  C„),  as  such 


/e[e(x,  y;  A,  b)  |iho,  Ax,  Ay]  =  ^e-|[t(x,,)-p(x,y|Ax,Ay)]cA[t(x,,)-p(x,,|Ax,A,)] 

where  C„  is  the  expected  sensor  noise  covariance  matrix  in  the  resulting  difference 
image.  Further,  let  {ax,  ay)  be  the  root  mean  squared  mis-registration  errors  and  as¬ 
sume  a  normal  prior  PDF  for  the  unknown  residual  spatial  mis-registration 


fAXAvi^^^^y) 


‘Zxcjx  (Ty 


— 


(15) 
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The  transformation  parameters  A  and  b  determined  by  chronochrome  (CC)  are 
assnmed  to  represent  the  maximnm  likelihood  estimates,  or  the  covariance  eqnal- 
ization  (CE)  estimates  are  nsed  as  a  snbstitnte.  The  change  test  statistic  can  be 
transformed  into  the  log-likelihood  by  d{x,y)  =  In  l{x,y).  Snbstitnting  the  results 


from  Eq.  (12)  to  Eq.  (15)  into  this  log-likelihood  test  statistic  and  ignoring  constant 


terms  that  are  not  data-dependent,  the  GLRT  test  statistic  becomes 


d{x,y)  =  min  \  [t{x,y)  -  p{x,y\Ax,  Ay)f  Cj[t{x,y)  -  p(x,  y|Ax,  Ay)]  +  ^  + 

A..A,  ^  y  j 

(16) 


Performing  the  minimization  dictated  by  the  GRLT  GD  statistic  in  Eq.  (16)  is 


straightforward  for  integer  pixel  mis-registration.  This  minimization  provides  the 
exact  same  result  as  the  GG  test  statistic  when  the  sample  covariance  matrix  of 
the  difference  image  Ce  is  used  as  an  estimate  of  the  noise  covariance  matrix  C„, 
but  with  a  minimization  over  different  test-to-reference  pixel  associations.  For  non¬ 
integer  pixel  mis-registration  increments,  some  interpolation  procedure  is  needed. 
In  the  general  case  where  the  mis-registration  could  be  more  than  one  pixel,  the 
minimization  needs  to  be  performed  over  integer  pixel  shifts  and  performed  over 
the  interpolation  between  integer  pixel  shifts.  The  assumption  is  that  the  image-to- 
image  registration  step  preceding  GD  can  provide  a  residual  mis-registration  better 
than  one  pixel  (which  is  typically  the  case)  and  only  explicitly  deal  with  the  case  of 
sub-pixel  shifts,  which  makes  the  mathematics  more  transparent.  It  should  be  noted, 
however,  that  the  method  developed  under  this  assumption  is  easily  generalized  to 
scenarios  with  mis-registration  errors  of  multiple  pixels  by  applying  it  over  multiple 
pixel  regions.  This  generalization  is  useful  for  extending  the  methods  described  to 
address  other  phenomenon  such  as  resampling  errors  introduced  during  image-to- 
image  registration  or  parallax  errors  caused  by  different  viewing  geometries  (addressed 
later  in  this  Ghapter). 
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The  prediction  p(a;,  ?/| Ax,  A?/)  is  estimated  in  two  dimensions  from  a  bilinear 
mixture  of  the  four  nearest  neighbors  of  p(x,|/).  The  Ax  and  Ay  provide  the  abso¬ 
lute  measures  of  the  mis-registration,  which  can  be  positive  or  negative,  and  include 
both  full  pixel  and  fractional  pixel  shifts.  However,  if  the  cost  function  is  limited  to 
fractional  mis-registration  and  adjusted  for  each  quadrant  around  the  pixels  then  a 
=  I  Ax  I  and  f3  =  |A|/|,  thus  the  estimate  is  given  by 

p(x,  y;  QlAx,  Ay)  =  (1  -  a)(l  -  /3)pJo  +  a(l  -  /5)p?o  +  (1  -  a)^P?i  +  a/5p?i-(17) 


where  Poo)P?0)Poi)  P?i  bilinear  mixtures  using  the  restricted  form  for 

quadrant  Q  (Table  around  a  center  pixel.  Note  in  Table  that  the  location  of  the 
four  pixels  mixed  in  the  bilinear  interpolation  are  merely  changed  to  maintain  a  and 
(3  between  zero  and  one. 

Table  3.  Definition  of  the  bilinear  mixture  spectra  in  each  of  the  four  quadrants. 


Quadrant  II  (Q=II) 

Poo  =  Ar(x,t/)  +  b 

Pio  =  Ar(x  -  1,1/)  -h  b 

Poi  =  Ar{x,y  +  1)  -h  b 

Pii  =  Ar(x  -  1,1/  -h  1)  -h  b 

Quadrant  I  (Q=I) 

Poo  =  Ar(x,t/)  -F  b 

Pio  =  Ar(x  -M,t/)  -h  b 

Poi  =  Ar{x,y  +  1)  -h  b 

Pii  =  Ar(x  -M,  ?/  -M)  -F  b 

Quadrant  III  (Q=III) 

Poo  =  Ar(x,t/)  -h  b 

Pio  =  Ar(x  -  1,1/)  -h  b 

Poi  =  Ar(x,//  -  1)  +  b 
pii  =  Ar(x  -  l,y  -  1)  +  h 

Quadrant  IV  (Q=IV) 

Poo  =  Ar{x,y)  -F  b 

Pio  =  Ar(x  -h  1,1/)  -h  b 

Poi  =  Ar(x,i/  -  1)  +  b 
pn  =  Ar(x  -h  1,1/  -  1)  b 

Dehning  q{x,  y;  Ax,  Ay)  as  the  function  within  the  braces  minimized  in  Eq. 


(16), 


q{x,  y,  Ax,  Ay)  =  [t(x,  y)  -  p{x,y\Ax,  Ay)f  ^  [t(x,  y) 


p{x,y\Ax,Ay)] 
Ax^  Aw^ 

H - n - 1 - 


(18) 


the  test  statistic  is  found  by  minimizing  q{x,y;a,  (3)  over  each  quadrant  Q=1...4 
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independently  subject  to  the  constraints,  0  <  a  <  1  and  0  <  /?  <  1. 

Inserting  Eq.  0  into  Eq.  (|I^,  then  multiplying  and  collecting  terms,  the  cost 
function  is  rewritten  as 


q{x,  y;  a,  /?,  Q)  =  +  b^a^jS  +  c^al3‘^  +  (i^a/3  + 


where  Q  again  denotes  the  quadrant  (note  that  the  /  above  is  a  coefficient  and  not  a 
PDF)  and  the  coefficients  of  the  cost  function  are  data-dependent  dehned  as 

=  [Pm  +  P?o-pS)-Pn]^C;^[p^i  +  pfo-pg)-p®] 

=  -2[p?o-Poo]^C;^[p^i  +  pfo-PTO-Pn] 

=  -2[p^i-p^o]^C;^[p^i  +  pfo-PTO-Pn] 

=  2[t{x,  y)  -  pS]'^C-1[pS  +  p?o  -  P?o  “  Pn]  +  2[p?o  “  P?o]^C-^[pJi  -  pg,] 
e'^  =  [p?o  -  [p%  -  p^o]  +  4 

^  X 

=  [pS-PTO]^C;^[pg-p^o]  +  ;4 

gQ  =  -2[t{x,  y)  -  Poo]^C;^[pfo  -  P^] 

=  -2[t{x,  y)  -  Poo]'^C;^[p^i  -  p^o] 

=  [t(x,|/)  -  pg)]^C;^[t(x,|/)  -  p^o] 


Since  Eq.  (19)  is  a  fourth  order  polynomial,  numerical  methods  are  required  to 
perform  the  constrained  minimization  in  each  quadrant.  Such  methods  depend  on 
computing  the  gradient  vector  and  Hessian  matrix,  dehned  as  (dropping  the  quadrant 
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index  Q  for  notational  ease) 


Vq  = 


da 

d/3 


(20) 


H 


d^q  d^q 
da^  dj3a 

d^q  d^q 
da/3  8/3'^ 


The  first  order  partial  derivatives  are  defined  as 


dq 

da 

dq 


2,(ia(3‘^  T  2,ba/3  T  c0^  dl3  T  260;  T  q  , 
2flo  /3  T  ba  T  2,ca/3  da  T  2,f  /3  h. 


The  second  order  partial  derivatives  are  dehned  as 


d‘^q 

da"^ 

d^q 

w 

d^q 
da  13 
d^q 
d(3a 


2af3^  +  2b/3  +  2e, 

2ao^  +  2co  +  2/, 

4(20/3  T  2ba  T  2c/3  T  ci, 
4(2o/3  T  2ba  T  2c/3  T  d. 


(21) 


(22) 


(23) 


The  minimization  problem  requires  a  suitable  optimization  algorithm  to  solve  but 
the  convergence  rate  can  be  slow;  hence  in  the  next  section  an  analytical  solution  is 
derived  to  support  an  inexpensive  computational  approximation. 
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3.1.2  Derivation  of  a  quadratic  approximation  to  the  GLRT  test  statis¬ 
tic 

A  second  order  Taylor  series  expansion  of  the  nonlinear  cost  function  q{x,  y,  a,  /3) 
can  be  used  to  support  an  analytical  minimization  and  closed-form  test  statistic  that 
should  support  a  more  computationally  efficient  change  detector.  For  each  of  the  four 
quadrants,  hrst  hnd  the  unconstrained  minimum  based  on  a  quadratic  approxima¬ 
tion  about  the  center  of  the  quadrant  (approximated  by  assuming  the  center  pixel), 
a  =  I3  =  ^.  If  this  unconstrained  minimum,  designated  as  {a*,  (3*),  falls  within  the 
constraint  boundaries  (0  <  a*  <  1  and  0  <  <  1),  it  is  accepted  as  the  constrained 

minimum.  Otherwise,  the  minimum  among  each  of  the  four  boundaries  of  the  quad¬ 
rant  is  determined,  and  the  smallest  value  of  the  four  is  accepted  as  the  constrained 
minimum.  No  approximation  is  required  along  these  quadrant  boundaries  as  the  cost 
function  is  quadratic  for  a  hxed  a  or  [3. 

The  two-dimensional  function  g(a;,  y]  a,  (3)  can  be  expanded  about  a  point  {d ^13') 
into  a  second  order  Taylor  series  as 


/  '  n'\  /  '  n'\ 

q{x,  y,a  ,/3)  =  q{x,  y,a  ,/3)+  — 


d^q 


dd^ 


1  /  'n2  , 


1  .  ,  <9? 

a  ,p 


/  nf  L\ 

a  ,p 


d^q 


dad  (3 


-{a-a){(3  -  (3)  +  ... 

'  a' 

CL  ,p 


(24) 

(25) 


Taking  the  derivatives  with  d  =  {3' 


and  keeping  only  the  terms  explicitly  shown 


46 


in  Eq.  (25),  the  quadratic  approximation  is 


q{x,y]  a,p') 


a  b  +  c  d  +  e  +  f  q  +  h 


+ 


a  +  c  b  +  d 
—r-  +  -^^  +  ^  +  9 


a--  1  + 


a  +  b  c  +  d 

—  +  —  +  /  +  /" 


a  b 
4  +  2+' 


a--  I  + 


a  c 

4  +  2+-'^ 


+  [a +  b  +  c  +  d][a  —  -ji(3  —  - 


(26) 


where  the  variables  a,  b,  c,  d,  e,  /,  g,  h,  k  are  dehned  in  Section  3.2.1.  Setting  the  deriva¬ 
tives  of  this  quadratic  approximation  with  respect  to  both  a  and  (3  to  zero,  the 


unconstrained  minimum  is  located  at  a*  and  f3*  per  Eq.  (27)  and  Eq.  (28). 


* _ {8d+8c+8b+8a)h-\-(—16f—8c—4:a)g+(4:C-\-8b+8a)f-\-(—4:C—2b—4:a)d—2c^-\-{—2b—3a)c—2b^—4:ab—2a^ 

{32e+l6b+8a)f+{l6c+8a)e—8d‘^  +  {—l6c—l6b—l6a)d—8c^  +  {—8b—l2a)c—8b'^  —  12ab—6a‘^ 


(27) 

(16e-|-8b-|-4a)ft-|-(— 8d— 8c— 8b— 8a)g-l-(— 8c— 4b— 8a)e-|-(2c-|-4&-|-4a)(i-|-2c^-|-(26-|-4a)c-|-2b^-|-3ab-|-2a^ 

(32e-|-16b-|-8a)/-|-(16c-|-8a)e— 8d2  +  (— 16c— 16b— 16a)ci— 8c2  +  (— 8b— 12a)c— 8b2  — 12ab— 6a2 


(28) 


The  minima  along  the  four  boundaries  are  computed  directly  per  Eq.  (23).  For 


(do,  0): 


For  (0,/3o): 


dq 

da 


n  .  y 

—  (J  — ^  dg  — - 


(oo,0) 


2e 


(29) 


dq 

da 


—  0  — )■  1^0  — 


(0,/3o) 


h 

2/ 


(30) 
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For  («!,  1): 


dq  ^  ^  (c  +  d  +  g) 

— —  —  (J  — y  Qfl  — - - ; - - 

(ai,i)  2{a  +  b  +  e) 


For  (l,/3i): 


^  =0  ^  B  jb  +  d  +  h) 

(1,^30  2(a  +  c  +  /) 


(31) 


(32) 


The  minimum  test  statistic  for  a  given  quadrant  Q  =  1 ...  4  is  then  given  by 


•Tim  {q{x,  y;  ag,  0,  Q),  q(x,  y;  0,  ^o,  Q),(l{x,  y;  Oi,  1,  Q),  q{x,  y;  1,  /3i,  Q)} 


0  <  a*  <  1  and  0  <  ^*  <  1 
otherwise 


(33) 


The  test  statistic  d{-)  is  then  the  overall  minimum  over  the  four  quadrants: 
d{x,y)  =  min  d{x,y;Q)  ioY  Q  =  1  ...  4. 

3.1.3  Comparison  with  local  co-registration  adjustment 


The  foundation  for  both  the  LCRA  and  the  GLRT  quadratic  approximation  is 
based  on  an  anomalous  change  detector.  The  LCRA  performs  an  interpolation  of 
the  test  statistic,  which  becomes  quadratic.  The  GLRT  results  in  a  fourth-order 
(quartic)  test  statistic  that  is  approximated  as  a  quadratic  using  a  second  order 
Taylor  series  expansion.  Both  algorithms  seek  to  improve  detection  performance 


for  mis-registered  data.  Unlike  LCRA  (Section  2. 4. 2. 3),  the  GLRT-based  approach 
developed  in  this  dissertation  incorporates  a  prior  probability  density  function  for  the 
residual  mis-registration,  which  might  be  critical  in  inhibiting  the  local  optimization 
from  hnding  minima  that  are  improbable  relative  to  the  statistics  of  the  actual  mis¬ 
registration  between  images.  The  GLRT  further  allows  one  to  incorporate  additional 
prior  information  regarding  the  image  mis-registration  statistics,  as  might  be  available 
to  compensate  for  other  image  acquisition  or  image  preprocessing  errors.  The  GLRT 
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formulation  results  in  a  closed  form  solution  for  the  test  statistic  that  is  minimized 
relative  to  the  unknown  local  shift  parameters.  One  major  downfall  of  LCRA  from  the 
perspective  of  this  research  is  that  it  does  not  provide  a  basis  to  extend  it  to  address 
the  parallax  problem.  The  GLRT  approach  provides  the  generality  to  incorporate  a 
priori  information  about  the  expected  direction  of  parallax  errors. 

3.2  Accounting  for  image  parallax 

To  account  for  parallax,  there  is  hrst  a  need  to  explicitly  extend  the  method  to 
handle  multiple  pixel  shifts.  In  this  case,  let 


Ax  =  m  +  a, 

A|/  =  h  +  /3. 


where  (m,  h)  are  whole  pixel  shifts  and  0  <  a  <  1  and  0  <  /?  <  1.  As  such,  the 


Quadrant  I  algorithm  is  used  in  all  cases  (see  Section  3.1.1)  with  the  addition  of  the 
(m,  h)  pixel  shift.  That  is,  the  sub-pixel  part  of  (Ax,  Ar/)  is  always  assumed  to  be  in 


the  positive  x  and  y  direction  from  the  full-pixel  {m^n)  part.  Following  Eq.  (10),  the 
prediction  then  becomes 


p(x,?/;  |h,m,Q;,/3)  =  (1  -  «)(!  -  /3)poo  +  «(1  -  /3)pio  +  (1  -  «)/3poi  +  a/3pii,  (34) 
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where,  per  Table 


Poo(a;,  y,  fh,  h)  =Ar(x  +  m,  ?/  +  h)  +  b, 

Pio(x,  I/,  m,  h)  =Ar(a;  +  m  +  1,  ?/  +  h)  +  b, 
Poi(a;,  y,  m,  h)  =Ar(a;  +  rh,  y  +  h  +  1)  +  b, 
Pii(a;,  y,  m,  h)  =Ar(a;  +  m  +  1,  y  +  h  +  1)  +  b. 


Note  that  fh  and  h  can  be  positive  or  negative. 

Let  0  be  the  direction  of  the  observed  parallax,  dp  the  RMS  parallax  error,  and 
am  the  RMS  residnal  mis-registration,  both  measnred  in  pixels.  The  expected  RMS 
parallax  error  can  be  estimated  by  the  change  in  viewing  geometry  and  the  3-D  scene 


strnctnre.  This  is  shown  in  Fig.  15  where  dr  is  the  occlusion  from  incidence  angle 
9r  and  dt  is  the  occlusion  from  incidence  angle  9t.  Thus  dt  —  dr  can  be  estimated 
as  the  parallax  error,  ap,  due  to  the  height,  h,  of  the  building  structure.  The  mis- 


Figure  15.  Image  depicting  geometry  of  acquiring  a  test  and  reference  image. 


registration  is  assumed  to  be  direction- independent  such  that  am  =  (^x  =  (^y  This 
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assumption  is  made  for  two  reasons.  First,  a  priori  estimates  of  the  misregistration 
parameters  ax  and  ay  are  unavailable  and  the  values  are  difficult  to  determine  without 
extensive  studies.  Second,  there  is  convenience  in  making  the  mathematics  cleaner 
given  the  assumption.  In  the  principal  point  coordinates  (A.^,  Arj)  where  A,^  is  along 
the  parallax  direction  and  Arj  is  orthogonal  to  it,  the  prior  PDF  for  the  local  mis¬ 
registration  is 


1  1 
(27r)  + 


Ari'^ 


2(ct^) 


(35) 


In  Fig.  [l6|(left),  the  Ax  and  Ay  axes  represents  the  zero  parallax  case  which  is 
rotated  about  the  origin  through  an  acute  angle  0  resulting  in  the  A^  and  Arj  axes 
corresponding  to  the  direction  of  parallax  (^)  and  the  orthogonal  direction  (p).  The 
angle  0  is  the  direction  of  parallax  relative  to  the  x-axis.  Thus,  a  given  point  p  has 
coordinates  (Ax,  Ay)  in  the  image  coordinate  system  and  (A,^,  Arj)  in  the  parallax- 
oriented  coordinate  system.  In  Fig.  [T^right),  the  line  from  the  origin  to  the  the 
point  p  has  magnitude  r  and  angle  6  in  the  (A^,  Ay)  coordinate  system.  Using  this 
observation,  the  relationship  between  A^  and  Ay  to  Ax  and  Ay  (as  illustrated  in 
Fig.  [l6|(right))  is  given  by 


A^  =  rcos6  Ay  =  rsm6, 

Ax  =  r  cos(0 -f  6*)  A?/ =  r  sin(0 6*). 


Using  the  addition  formula  for  the  cosine  and  sine  functions.  Ax  and  Ay  are  expressed 
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Figure  16.  The  relationship  between  coordinates  as  they  describe  the  same  point,  p  in 
the  plane,  (left)  The  Ax  and  Ay  axes  have  been  rotated  about  the  origin  through  an 
angle  (p  (the  parallax  direction  angle)  to  produce  the  Af  and  A77  arxes.  (right)  A^  and 
Arj  are  related  to  Ax  and  Ay  through  the  blue  line  from  the  origin  to  point  P. 


in  terms  of  and  Arj  as 


Ax  =  r  cos(0  +  9)  =  r(cos0cos6*  —  sin0sin6*), 

Ax  =  (r  cos  9)  cos  0  —  (r  sin  9)  sin  0, 

Ao:  =  cos  0  —  At]  sin  0.  (36) 


and 


Ay  =  r  sin(0  +  9)  =  r(sin0cos6*  +  cos  0  sin  0), 

Ay  =  (r  cos  9)  sin  0  —  (r  sin  9)  cos  0, 

At/ =  A^  sin  0  +  At]  cos  0.  (37) 


Solving  for  A^  and  Ay  in  Eq.  (36)  and  Eq.  (37)  yields 


A,^  =  Ax  cos  0  +  Ay  sin  0, 
Af]  =  —Ax  sin  0  +  Ay  cos  0. 


(38) 

(39) 


After  expressing  the  new  coordinates  from  axis  rotation,  an  evolution  of  the  normal 
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prior  PDF  is  formulated  to  include  an  expected  parallax  orientation  (0)  and  RMS 


parallax  error  (dp).  Thus,  substituting  Eq.  (38)  and  Eq.  (39)  into  Eq.  (35)  yields 


fAxAvi^^y^y)  = 


(Ax  COB  tl>+ Ay  Bin  cji)^  _  (- Aa:  sin  <^+At/ cos  </))2 


(27r)  + 


2(ct4i  +  o'p) 


m'^pJ  g 


(40) 


The  original  cost  function  in  Eq.  ( 18 )  then  becomes 


q{m,n,a,/3)  =  [t{x,y)  -  p{x ,  y\Ax ,  Ay)Y C ^\t{x ,  y)  -  p{x,y\Ax,  Ay)] 

I  (Ax  COS  Ay  sin  0)^  ,  (— Axsin0+Aycos0)^ 


(41) 


For  each  (m,  n)  the  polynomial  is  of  the  same  form  as  before  but  the  penalty  terms 
are  different.  They  are  found  by  expanding  two  quadratic  equations 


{Ax  cos  0  +  Ay  sin  0)^, 

(42) 

—Ax  sin  0  +  Ay  cos  0)^. 

(43) 
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The  coefficients  used  in  the  cost  function  described  in  Eq.  ( 19 )  now  become 


a  = 


b  = 


c 

d 


[Poi  +  PlO  ~  Poo  ~  Pll]^C„^[poi  +  Pio  —  Poo  ~  Pll] 
-2[pio  -  Poo]^C“^[poi  +  Pio  -  Poo  -  Pll] 

-2[poi  -  Poo]^C“^[poi  +  pio  -  Poo  -  Pi 


ni 


=  2[t{x,y)  -  Poo]^C„^[poi  +  pio  -  Poo  -  Pii]  +  2[pio  -  Pooj"  C„^[poi  -  Poo] 


Tr^-u 


cip  2  cos  0  sin  ( 

(T^  (T^  +  CT^ 


f 

9 


—  [Pio  —  Poo]'^C„^[pio  —  Poo]  + 
=  [Poi  -  Poo]^C~^[poi  -  Poo]  + 


+  ClpSin^0 

+  ^p) 

+  tX^COsV 

^  m\^  m  '  ^pJ 


=  -2[t(a;,i/) -poo]^C„^[pio-poo] 

2  cos  0(m  cos  (j)  +  n  sin  0)  ^  2  sin  0(m  sin  (p  —  n  cos  0) 


h  = 


2 
P 

Tr^-lr 


at 


-2[t(x,?/)  -  Poo]  [poi  -  Poo] 

2  sin  0(m  cos  0  +  h  sin  0)  2  cos  0(m  sin  cp  —  h  cos  (p) 


k  = 


0-2  _j_  02 

^  m  '  ^  p 


(7 


.  ,  ,  vJv-.-ir  /  N  1  (mcos0  +  nsin^)^  (msin^  —  hcos^)^ 

[t(x,  y)  -  Poo]  C„  [t(x,  y)  -  poo]  + - ^ ^ - + 


at 


Calculating  the  parallax  direction  through  stereo  correspondence  was  described 
previously  in  Section  2.5[  The  stereo  correspondence  method  utilized  in  this  work  is 
based  on  matching  SIFT  descriptors.  The  SIFT  descriptors  are  computed  and  com¬ 
pared  using  a  similarity  measure  to  determine  correspondence  points.  The  similarity 
is  dehned  as  the  cosine  of  the  angle  between  SIFT  vectors  in  the  reference  image  and 
test  image  as 


sim{Yir,Vt)  = 


Pr  Pt 


(44) 


iPr-ll  IlPdl 

where  p^  and  p*  are  SIFT  descriptor  vectors  from  the  reference  and  test  images,  re¬ 
spectively,  and  1 1  ■  1 1  denotes  the  Euclidean  norm.  Matching  descriptor  points  between 
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the  reference  and  test  images  is  accomplished  by  computing  the  similarity  measure 
between  all  feature  pairs  where  the  two  feature  pairs  are  declared  a  match  when  their 
similarity  is  within  a  threshold.  For  a  set  of  correspondence  points,  the  maximum 
is  the  most  likely  shift  for  parallax.  This  eliminates  small  shifts  likely  from  slight 
mis-registration.  The  parallax  direction  is  computed  in  terms  of  the  coordinate  shift, 
expressed  in  pixels,  by  selecting  the  correspondence  points  with  the  largest  Euclidean 
distance  between  index  values, 

k  =  argmax  {xt^i  -  Xr,if  +  {yt,i  -  yr,if  ,  (45) 

i 

where  i  is  the  pair  of  correspondence  points  under  consideration,  {xr^i,yr^i)  and  {xt^i,  yt^i) 
are  the  indices  of  correspondence  point  i  in  the  reference  and  test  images  respectively, 
and  k  is  the  reference  and  test  correspondence  point  used  to  determine  the  parallax 
angle.  Next,  the  parallax  direction  angle  (0)  (shown  in  Fig.  [T^  for  the  given  corre- 


Figure  17.  Calculating  the  parallax  direction  angle,  cj),  from  the  calculated  correspon¬ 
dence  point  shift. 

spondence  point  is  calculated  as 

i  ,  (  yt,k  yr,k  \ 

0  =  atan  -  . 

\^t,k  ^r,k  J 
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If  the  denominator  is  zero  then  0  =  90°. 


The  correspondence  points  allow  a  parallax  direction  angle  estimate,  0,  then  to 


be  used  in  Eq.  (41)  for  computation  of  the  test  statistic. 


3.3  Summary 

This  chapter  dehned  two  algorithms  used  to  address  the  problems  of  image  mis¬ 
registration  and  image  parallax  as  they  pertain  to  the  CD  problem.  The  mis-registration 
algorithm  incorporates  a  spatial  mis-registration  model  into  a  GLRT-based  change 
detector  leading  to  a  fourth  order  polynomial  test  statistic.  A  second  order  Taylor 
series  expansion  of  the  nonlinear  cost  function  is  then  used  to  derive  a  closed-form 
test  statistic  to  support  a  more  computationally  efficient  change  detector.  The  second 
algorithm  leveraged  the  inherent  relationship  between  the  reference  and  test  image 
views  and  pixel  correspondence  matching  to  perform  parallax  mitigation.  This  re¬ 
quires  a  search  for  correspondence  points  in  order  to  calculate  the  parallax  angle,  0, 
to  determine  if  in  the  pixel-level  change  detector,  a  true  change  pixel  exists  even  along 
the  parallax  direction. 

Chapter  |IV|  describes  the  test  data  used  in  algorithm  analysis.  Algorithm  perfor¬ 
mance  of  both  proposed  algorithms  are  compared  to  existing  anomalous  CD  schemes. 
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IV.  Experimentation  and  Results 


4.1  Overview 


The  previous  chapters  presented  necessary  background  material  and  derived  novel 
hyperspectral  change  detection  (HSCD)  algorithms.  In  this  chapter,  hyperspectral 
data  used  to  detail  the  algorithms  are  presented  and  described.  The  generalized 
likelihood  ratio  test  (GLRT)  mis-registration  compensation  algorithm  presented  in 


Chapter  III  is  compared  with  current  state-of-the-art  anomalous  change  detection 
(CD)  algorithms,  namely  chronochrome  (CC)  and  covariance  equalization  (CE),  as 
well  as  the  local  co-registration  adjustment  (LCRA).  The  GLRT-based  parallax- 
compensation  algorithm  is  compared  to  CC,  and  parameter  sensitivity  is  explored 
for  both  GLRT-based  algorithms.  All  of  the  results  presented  in  this  chapter  were 
developed  and  applied  in  the  MATLAB®  software  environment. 


4.2  Data  description 

There  are  two  data  sets  used  to  evaluate  the  mis-regristration  algorithm,  the  Air 
Force  Research  Lab  (AFRL)  tower  data  and  the  Civil  Air  Patrol  (CAP)  Airborne 
Real-time  Cueing  Hyperspectral  Enhanced  Reconnaissance  (ARCHER)  [76]  Mojave 
data  sets.  In  particular,  a  simulated  change  to  the  CAP  ARCHER  eg  Mojave  data 
set  is  used  to  compare  mis-registration  compensating  change  detector  results  with 
LCRA.  An  additional  data  set  is  presented  and  is  used  for  evaluating  the  parallax 
compensation  algorithm,  a  synthetic  data  set  generated  by  the  Digital  Imaging  and 
Remote  Sensing  Image  Generation  (DIRSIG)  system  [77]. 
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4.2.1  Data  to  evaluate  mis-registration-compensating  change  detec¬ 


tion  algorithm 


4.2. 1.1  Tower  data 


The  tower  data  shown  in  Fig.  used  for  testing  the  proposed  detector,  was 
collected  in  a  controlled  manner  from  a  tower  using  a  visible  to  near  infrared  (VNIR) 
hyperspectral  imaging  spectrometer  mounted  on  a  pan  and  tilt  assembly.  The  spectral 
range  of  the  imager  is  0.46  —  0.9  /im  and  the  spatial  resolution  is  approximately  4  cm. 
The  tower  data  was  collected  at  Wright  Patterson  Air  Force  Base,  Ohio.  While  the 
data  collection  effort  extended  from  August  2005  through  May  2006,  this  particular 
pair  was  collected  on  October  14,  2005  near  solar  noon  [TS].  Less  than  one  half 
hour  passes  between  collecting  the  test  and  reference  image  pair  thereby  limiting 
illumination  and  vegetation  changes  between  images.  In  doing  so,  any  degradation 
in  performance  that  may  occur  can  be  attributed  mostly  to  prediction  error  due  to 
artihcially  introduced  mis-registration  rather  than  errors  introduced  by  shadowing  or 


vegetation  changes  d-  The  top  two  images  in  Fig.  ^  are  a  color  composite  of 
the  data.  The  scene  is  composed  of  mostly  grass  and  trees  with  coated  aluminum 
panels  placed  along  the  treeline.  One  change  of  interest  is  identihed  in  the  scene, 
corresponding  to  a  small  green  tarp  bundle  added  to  the  test  image.  The  challenge 
for  the  CD  algorithm  is  to  detect  the  small  change  in  the  presence  of  mis-registration. 
A  pixel-level  truth  map  for  the  spatial  pixel  dimensions  (126x126)  tower  data  is 
extracted  manually  based  on  the  CC  test  statistic  to  evaluate  the  results  and  is 


displayed  in  Fig.  18  (bottom).  The  truth  image  represents  pixels  as  black  (0)  for 


no  change  and  white  (1)  for  true  change.  The  change  target  in  the  truth  mask  is 
a  rectangular  grid  of  3x5  pixels  with  a  one  pixel  ignore  mask  along  the  outer  edge 
of  the  change  target.  An  ignore  mask  identihes  pixels  around  the  periphery  of  the 
change  target  that  are  ignored  during  scoring  to  keep  them  from  biasing  the  results 
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because  the  ground  truth  for  these  pixels  is  uncertain.  The  background  consists  of 
15,841  pixels. 


Figure  18.  Change  pair  tower  data  set  used  for  mis-registration  analysis  (top).  These 
false  color  composites  are  produced  from  the  hyperspectral  data  cubes  by  selecting 
respective  red,  green  and  blue  bands.  Binary  target  truth  mask  for  the  tower  data  set 
(bottom). 

To  test  the  proposed  algorithms,  the  test  image  is  artihcially  mis-registered  by 
translating  the  original  high  resolution  data  then  mean  hltering  to  simulate  fractional 
pixel  shifts  as  shown  in  Fig.  [T^  The  test  data  is  shifted  one  pixel  in  both  the  positive 
X  and  negative  y  direction  with  respect  to  the  reference  spatial  scene  to  simulate  mis¬ 
registration.  A  sub-pixel  mis-registration  was  introduced  by  resampling  both  images 
(reference  and  test)  with  a  2x2  pixel  mean  hlter  to  create  the  sub-pixel  mis-registration 
shift.  In  order  to  eliminate  edge  artifacts,  the  images  are  cropped  to  ensure  the  same 
number  of  samples  is  used  to  estimate  test  statistics.  Shifting  integer  pixel  amounts 
for  the  original  high  resolution  imagery  and  then  synthetically  reducing  resolution  by 
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averaging  full  pixel  blocks  avoids  potential  spectral  interpolation  artifacts. 


(1+244+6)/4 
=  3.25 

4.75 

(3+840+0)/4 
=  2.75 

5.5 

7.5 

4.25 

3 

4 

2.25 

Figure  19.  An  example  of  the  simulated  mis-registration  process.  The  images  on 
the  left  represent  high  resolution  imagery.  The  bottom-left  image  is  shifted  prior  to 
performing  block  averaging.  The  images  on  the  right  represent  the  new  images  at  a 
coarser  spatial  resolution  where  the  bottom  image  also  reflects  a  simulated  1/2  pixel 
shift  in  the  positive  x  and  negative  y  directions. 


4. 2. 1.2  Mojave  data 


CAP  ARCHER  [76]  hyperspectral  images  were  captured  at  the  collection  con¬ 
ducted  at  the  National  Test  Pilot  School  in  Mojave,  California.  The  objective  of  the 
collection  was  to  demonstrate  the  operational  utility  and  baseline  the  performance 
of  a  VNIR  hyperspectral  sensor  for  a  search  and  rescue  mission.  The  images  with  a 


spatial  pixel  dimension  size  of  238x208  are  displayed  in  Fig.  ^  Again,  a  synthetic 
mis-registration  is  introduced  in  order  to  perform  an  analysis  on  the  mis-registration 
compensation  algorithm.  The  same  approach  for  synthetically  mis-registering  the 
data  is  used  for  this  change  pair  as  was  described  for  the  tower  data.  Specihc  spatial 
portions  of  the  synthetically  mis-registered  data  were  swapped  to  simulated  change 
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occurring  between  the  reference  and  test  images.  Fig.  20  shows  a  synthetic  change 
for  testing  the  algorithms  that  is  created  by  interchanging  a  square  9x9  block  of  pix¬ 
els  from  one  spatial  location  in  the  image  to  another  (see  red  circled  image  chips  in 
Fig.  [fright).  This  provides  known  true  change  pixels  for  algorithm  evaluation  which 
can  be  seen  in  the  truth  mask  (bottom)  where  target  pixels  are  represented  as  white. 


Figure  20.  Top  images  show  CAP  ARCHER  data  set  with  synthetic  change.  The 
swapped  pixels  are  shown  in  red  circles  on  the  right  image.  Bottom  image  shows  the 
target  truth  mask. 


4.2.2  Data  to  evaluate  parallax-compensating  change  detection  algo¬ 
rithm 

4. 2. 2.1  Synthetic  DIRSIG  data 

The  parallax  GLRT  algorithm  is  tested  on  synthetic  imagery  generated  using  the 
DIRSIG  system  [77].  DIRSIG  is  a  software  suite  started  in  the  late  1980s  at  the 
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Rochester  Institute  of  Technology  (RIT)  in  Rochester,  NY.  The  program  predicts 
imager  radiance  values  using  a  physics  model  and  ray  tracing.  The  physics  model 
takes  into  account  thermal  characteristics,  atmospheric  properties  and  distortion,  HS 
imager  features  and  noise,  the  sensor  platform,  and  the  surface  composition  m 
The  DIRSIG  modeler  uses  actual  end-member  spectra  for  in-scene  materials,  and 
calculates  the  spectral  response  as  seen  by  the  sensor  for  the  specihc  environmental 
conditions.  The  software  is  used  by  civilian,  academic,  and  defense  agencies  for 
various  scientihc  and  engineering  studies  [2l  12^  1801-183] .  The  use  of  synthetic  imagery 
provides  the  advantage  of  having  complete  truth  information  and  exact  positions  of 
the  change  targets  within  the  scene. 

An  artihcial  scene  that  mimics  an  Airborne  Visible  InfraRed  Imaging  Spectrom¬ 
eter  (AVIRIS)  sensor  with  224  spectral  channels  ranging  from  approximately 
400  -  2500  nm  is  constructed  to  represent  an  urban  scene  with  elevated  structures. 
The  spatial  resolution  is  approximately  0.3m  and  the  image  spatial  size  is  311x259. 
The  synthetic  scene  was  constructed  to  introduce  image  parallax.  The  rendered  scenes 
were  simulated  to  collect  at  10,000  ft  with  two  different  viewing  angles:  10°  and  20°. 
The  reference  and  test  images  are  formed  with  a  vehicle  departure  change  occurring 


between  instances  (red  circles,  top  of  Fig.  21). 

The  parallax  direction  is  estimated  by  manually  evaluating  the  building  edges. 
The  different  off-nadir  view  angles  between  the  reference  and  test  images  results  in 
parallax  and  this  is  the  manually  estimated  parallax  direction  from  “building  lean” 
where  the  sides  of  the  buildings  become  more  visible  with  an  increase  in  off-nadir 
viewing  angle.  The  manually  estimated  parallax  direction  average  is  approximately 


4  pixels  in  the  upward  direction  (90°,  see  Fig.  22). 

In  addition,  the  synthetic  data  allows  a  controlled  environment  for  knowable 
ground  truth.  DIRSIG  supports  per-pixel  truth  maps  that  provide  important  im- 
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age  formation  parameters  |77|.  The  true  parallax  in  this  imagery  can  be  calculated 
from  the  per-pixel  truth  maps  provided  by  DIRSIG.  The  displacement  or  true  pixel 
shift  (tps),  in  pixels,  can  be  determined  by  knowing  the  collection  angle,  the  height 


of  the  buildings,  and  the  ground  sampling  distance  (GSD)  (Fig.  15).  Therefore,  a 


truth  mask  is  developed  from  the  pixel  truth  information  furnished  by  DIRSIG  by 


calculating  Eq.  (47)  at  every  pixel 


tps  = 


h  tan  Or  h  tan  Of 


GSD 


GSD  ’ 


(47) 


where  h  is  the  building  height  and  Or  and  Of  are  the  reference  and  test  viewing  angles. 


The  parallax  truth  mask  is  shown  in  Fig.  23,  where  the  colorbar  on  the  right 
reflects  the  true  value  of  the  local  shift.  There  are  a  total  of  1,562  true  parallax 
pixels.  An  ignore  mask  was  created  for  pixels  above  the  ground  plane  but  not  on 


the  edge  of  a  structure  and  are  the  gray  pixels  in  Fig.  .  The  buildings  where 
the  parallax  does  not  occur  are  considered  to  be  ignore  pixels  because  although  the 
building  reveals  a  shift  in  the  parallax  truth  map,  the  algorithm  sees  ‘no  change’  as  it 
is  still  of  the  same  material.  The  ignore  pixels  account  for  8,807  pixels  out  of  the  total 
80,549  pixels  in  the  image.  The  background  pixels  (ground  plane  pixels)  experience 
no  height  difference  thus  is  zero  and  account  for  70,180  pixels. 
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Figure  21.  Synthetic  reference  image  of  an  urban  scene  with  elevated  structures  and  a 
stationary  red  vehicle  (red  circle)  taken  at  a  10°  viewing  angle  is  displayed  upper  left. 
On  the  upper  right  is  the  synthetic  test  image,  captured  from  an  alternate  angle  (20° 
viewing  angle)  with  the  vehicle  no  longer  in  scene. 
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Figure  22.  In-scene  building  edge  from  the  synthetic  reference  image  (collected  at  a 
10°  viewing  angle)  used  to  manually  estimate  the  parallax  direction  between  images 
(left  images).  Synthetic  test  image  (right  images)  of  the  same  in-scene  building  but 
captured  from  an  alternate  viewing  angle  of  20°  displaying  an  increase  in  visibility  of 
the  building  side  with  an  appraised  parallax  shift  of  4  pixels  in  the  upward  direction 
(90°).  Labels  on  both  image  chips  are  pixel  values. 


4.3  Algorithm  implementation  details 

4.3.1  Mis-registration  algorithm  implementation 

Both  the  CC  and  CE  anomalous  change  detectors  rely  on  a  linear  prediction  of 
the  reference  image  using  the  statistics  from  the  test  image.  Assuming  this  linear 
relationship  exists,  a  gain  matrix  (A)  and  offset  vector  (b)  are  estimated  and  applied 
to  the  reference  image  to  match  conditions  of  the  test  image.  Upon  applying  the 
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Figure  23.  Depiction  of  parallax  truth  mask.  The  parallax  effects  reside  along  the 
building  edges.  The  gray  pixels  reflect  ignore  masks  so  that  top  pixels  of  elevated 
structures  are  ignored  during  analysis. 


predictor,  the  difference  image  between  the  prediction  and  test  images  is  obtained, 
followed  by  an  anomaly  detector  to  determine  if  a  change  has  occurred.  This  is  the 


typical  process  for  anomalous  CD  and  is  shown  in  the  dotted  box  in  Fig.  24  The 


remaining  blocks  in  Fig.  24  show  the  CD  process  based  on  the  derived  GLRT  in 


Section  3.1  where  the  minimization  block  corresponds  to  minimizing  the  objective 


function  (described  in  Eq.  (19)  through  a  numerical  optimization  technique  or  the 
quadratic  approximation  in  Eq.  ([26|).  This  minimization  is  performed  for  each  of  the 
four  quadrants  (refer  to  Table  where  the  minimum  over  the  four  quadrants  is  the 
minimized  test  statistic.  In  Fig.[^  C„  is  the  covariance  matrix  of  the  difference  image 
based  on  the  CC  predictor  and  ax,  Cy  are  user-specified  parameters  characterizing  the 
root  mean  squared  (RMS)  mis-registration. 
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Detection 
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Figure  24.  Flow  chart  for  applying  the  mis-registration  GLRT  change  detection  algorithm. 


The  minimization  problem  defined  in  this  dissertation  requires  a  suitable  opti¬ 
mization  algorithm  to  solve.  In  general,  there  are  two  categories  of  approach:  deriva¬ 
tive  and  non  derivative.  Derivative  approaches  are  based  on  either  the  gradient  or 
Hessian  (first  and  second  derivatives  respectively).  Gradient-based  approaches,  such 
as  steepest /gradient  descent  (a  greedy  method),  are  popular  due  to  ease  of  prob¬ 
lem  formulation  but  tend  to  converge  slowly  since  movement  towards  the  (possi¬ 
bly  local)  minimum  is  linear.  Convergence  rate  is  improved  through  various  algo¬ 
rithmic  enhancements,  such  as  the  addition  of  a  momentum  term  in  back  propa¬ 
gation  neural  networks  [HI].  Hessian-based  approaches  generally  follow  Newton’s 
method  (see  e.g.^  [HS]),  which  has  a  quadratic  convergence  rate  and  hence  converges 
much  quicker  than  gradient-based  approaches.  Variants  of  Newton’s  method  {e.g.^ 
Broyden-Fletcher-Goldfarb-Shanno  [86|  and  Gauss-Newton  IH71)  are  used  to  overcome 
certain  limitations  such  as  the  requirement  of  the  Hessian  being  positive  definite  and 
non-singular.  The  nature  of  the  problem  may  require  one  to  seek  non  derivative- 
based  approaches.  Such  search  algorithms  include  the  class  of  genetic  algorithms  and 
stochastic  algorithms  {e.g.,  Simulated  Annealing  |HH])  to  name  a  few.  The  search 
literature  is  vast  and  myriad  good  options  exist,  each  with  its  merits  and  demerits. 
A  constrained  variant  [89]  of  the  Nelder-Mead  simplex  method  [90]  (a  direct  search 
approach)  that  transforms  the  constraints  and  uses  the  standard  unconstrained  form 
for  the  minimization  task,  is  used  in  this  dissertation. 

For  a  function  with  n  =  2  variables,  the  Nelder-Mead  method  creates  a  simplex 
of  n  -|-  1  points  (a  triangle  for  the  spatial  problems  explored  in  this  research),  and 
performs  a  search  that  compares  function  values  at  each  of  the  vertices.  The  vertex 
where  the  function  takes  its  largest  value  is  rejected  and  replaced  with  a  new  vertex.  A 
new  triangle  is  formed  and  the  search  is  continued.  The  process  generates  a  sequence 
of  triangles,  for  which  the  function  values  at  the  vertices  decrease  and  the  size  of  the 


triangles  is  reduced  locating  the  coordinates  of  the  minimum  point  pT], 


4.3.2  Parallax-compensation  algorithm  implementation 


The  parallax-compensation  flow  diagram  in  Fig.  25  is  similar  to  that  of  the  mis¬ 


registration  flow  diagram  in  Fig.  in  Section  [4 . 3 . 1 1  where  the  anomalous  CD  process 
remains  unchanged  but  the  minimization  process  is  noticeably  different.  In  Fig.  2^ 
the  minimization  block  minimizes  the  objective  function  with  the  parallax-derived 


coefficients  described  in  Section  |3.2|  to  include  integer  pixels  shifts. 

This  minimization  is  performed  for  multiple  (m,  n)  values  where  (m,  n)  is  estab¬ 
lished  based  on  (j).  The  search  region  was  set  manually  for  this  experiment  based  on 
the  fact  that  (j)  is  in  the  vertical  direction.  Thus,  rh  is  small  and  restricted  to  (—1,  0) 
while  h  is  bounded  by  an  assumed  <Jp  =  4  pixels  (— 2cTp  <n<  2ap). 

Note  that  the  search  region  can  very  easily  be  established  over  a  square  search 


region  bounded  by  —2ap  <  {fh^h)  <  2ap  (illustrated  in  Fig.  26),  however  the  prior 
PDF  allows  the  cost  function  with  lower  cost  in  the  0  direction.  Therefore,  choosing 
the  search  region  according  to  the  (j)  estimate  would  be  appropriate  (see  Fig.  26 
illustration  for  ‘low-cost’  and  ‘high-cost’  search  direction). 
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Figure  25.  Flow  chart  for  applying  the  parallax-compensation  GLRT  change  detection  algorithm. 


Test  Image 


Figure  26.  Parallax-compensation  algorithm  search  space  defined  for  parameters  (m,n). 


4.4  Scoring  methodology 


Pixel-level  scoring  (described  in  Section  2. 4. 5. 2)  is  used  to  generate  the  ROC 


curves  describing  pixel-level  success.  For  a  CD  system,  the  ROC  analysis  consists  of 
measuring  the  binary  response  of  the  detection  system  by  calculating  the  probability 
of  detection  as 


Pd  =  (48) 

where  D  are  properly  matched  pixels  and  T  are  true  change  target  pixels  and  the 
probability  of  false  alarm  with 


PpA  — 


M  -  D 


N-T 


(49) 


where  N  is  the  total  number  of  pixels  in  the  image  and  M  is  the  declared  changed  pix¬ 
els  at  a  particular  threshold.  The  {PfAi  Pd)  pah  is  one  point  on  the  ROC  curve.  This 
process  is  repeated  for  multiple  thresholds  where  M  and  D  will  vary  with  different 
thresholds. 
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At  each  decision  threshold,  and  PpA  are  calculated  by  comparing  a  binary 
detection  image  with  a  ground-truth  mask.  In  practice,  the  simplest  way  to  compare 
the  detection  image  and  the  truth  mask  is  to  make  a  pixel-level  comparison  exactly 
htting  the  dehnitions  of  Pp  and  Pfa-  The  manually  derived  ground-truth  mask  can  be 
difficult  to  dehne  and  skew  performance  if  not  accurate.  In  particular,  edge  pixels  may 
be  either  uncertain  or  contain  a  mixture  of  target  and  background,  especially  if  mis¬ 
registration  errors  are  present.  Thus,  after  deriving  the  ground-truth  pixels  based  on 
the  baseline  CC  test  statistic,  the  change  target (s)  are  surrounded  with  non-penalizing 
pixels  around  the  change  target  borders  (the  ignore  mask).  The  performance  is  further 
examined  when  the  ignore  mask  extends  into  the  target  mask  to  reduce  the  effect  of 
uncertain  edge  pixels. 

4.5  Experiments  and  results 

Performance  analysis  occurs  in  two  ways  in  this  dissertation.  First  is  a  visual 
analysis  of  the  detection  image.  The  visual  analysis,  when  combined  with  the  truth 
mask  and  knowledge  of  the  CD  scenario,  helps  provide  insights  into  what  the  change 
detector  is  doing  on  the  given  set  of  imagery.  It  further  provides  visual  conhrmation 
that  the  algorithms  are  performing  in  a  manner  that  makes  logical  sense.  Second 
is  a  quantitative  analysis  by  way  of  the  ROC  curve  and  the  AUC  derived  from  the 
ROC  curve.  The  ROC  curve  is  the  standard  method  of  performance  assessment  in 
the  detection  literature  and  is  used  in  this  dissertation.  In  the  case  there  are  non¬ 
dominating  ROC  curves,  the  AUC  provides  another  measure  to  comapre  performance. 
An  exception  is  for  the  DIRSIG  imagery,  where  the  Signal-to-Clutter  Ratio  (SCR)  is 
used  to  show  how  the  clutter  from  the  edge  artifacts  is  reduced  by  increasing 

Other  performance  analysis  approaches  exist  and  include  histogram-based  meth¬ 
ods  and  simply  plotting  intensity  versus  a  1-D  pixel  index.  A  common  histogram- 
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based  method  uses  a  histogram  of  the  change  image  to  determine  if  modes  exist  within 
the  detection  image  that  can  be  post-processed  to  identify  change  versus  no  change. 
The  latter  approach  is  not  as  common  and  may  be  difficult  to  interpret.  These  are 
not  pursued  in  this  dissertation. 


4.5.1  Mis-registration  algorithm 


This  section  describes  the  results  from  the  mis-registration  GLRT  algorithm  de¬ 


rived  in  Section  3.1.1  almost  exclusively  on  the  tower  data  before  considering  the 
Mojave  data.  First,  using  the  CC  predictor,  the  nonlinear  optimization  results  are 
analyzed  to  provide  evidence  that  the  algorithm  and  its  quadratic  approximation  are 
functioning  as  derived.  Second,  using  the  CC  predictor,  the  quadratic  approximation 
results  are  compared  to  the  nonlinear  optimization  approach  to  demonstrate  the  ef- 
hcacy  of  the  estimate.  Third,  using  the  CE  predictor,  performance  of  the  quadratic 
approximation  is  compared  to  the  CC  results.  Fourth,  the  quadratic  approxima¬ 
tion  using  the  CC  predictor  is  compared  with  LCRA  [33]  •  Finally,  model  parameter 
performance  impact  is  explored  using  the  quadratic  approximation. 


4.5. 1.1  Nonlinear  optimization 


An  inspection  of  the  unthresholded  detection  image,  using  the  CC  predictor, 
for  the  mis-registration  CLRT  algorithm  using  nonlinear  optimization  is  shown  in 


Fig.  27  After  experimenting,  the  RMS  mis-registration  parameters  were  chosen  as 
(Tx  =  CTy  =  0.05  for  computing  the  test  statistic  image  results  for  the  tower  data  set. 
The  arrows  indicate  the  mis-registration  artifacts;  the  test  statistic  is  visibly  lower  in 
Quadrant  II  which  is  expected  because  the  mis-registration  shift  is  in  that  direction. 
A  common  scale  is  used  where  black  and  white  pixels  in  the  change  mask  correspond 
to  a  test  statistic  dxnin{Qi,Q2,Q3,Q4,}  2/)  56.9  and  djjiax{Qi,Q2,Q3,Q4,}  2/)  745.5, 
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respectively.  The  detection  image  can  take  values  in  the  range  [0,oo).  Experiments 
show  that  Quadrant  II  suppresses  the  edge  effects  due  to  the  mis-registration  and 
is  the  correct  direction  in  which  the  actual  mis-registration  shift  of  Ax  =  1/2  and 
Ay  =  —1/2  occurred.  Therefore,  by  choosing  the  minimum  test  statistic  over  the 
four  quadrants  detected  change  caused  by  mis-registration  is  reduced. 

The  uncompensated  CC  test  statistic  image  (bottom  left),  test  statistic  using 
minimum  nonlinear  optimization  (top),  and  the  test  statistic  using  the  known  mis¬ 


registration  (bottom  right)  are  compared  in  Fig.  28  Uncompensated  CC  (Fig.  28 


bottom  left)  highlights  the  impact  of  the  mis-registration  effect  on  CD.  This  impact 
is  indicated  by  the  two  bright  white  vertical  lines  on  the  left  and  right  side  of  the  panel 
that  result  in  a  false  declaration  of  change.  When  using  the  GRLT  with  nonlinear 
optimization  (Fig.  |^top),  there  is  visual  evidence  that  false  changes  caused  by  mis¬ 
registration  are  alleviated  (note  the  darker  appearance  corresponding  to  a  lower  test 
statistic,  indicative  of  being  non  change  pixels).  The  test  statistic  using  the  known 


mis-registration  (Fig.  28  bottom  right)  represents  the  minimum  of  the  average  pixel 
spectra  for  each  of  the  four  quadrants  as  input  to  the  CC  algorithm.  The  edges 
here  offer  a  lower  test  statistic  and  are  accomplished  using  the  known  (true)  mis¬ 
registration. 


ROC  curves  are  shown  in  Fig.  29  and  it  is  evident  that  the  false  alarm  rate  is 
reduced  using  the  proposed  GLRT-based  method.  For  similar  false  alarm  rates,  the 
GLRT-based  approach  provides  two  to  three  times  the  detection  performance.  For 
similar  detection  results,  the  GLRT-based  approach  provides  one  to  two  orders  of 


magnitude  reduction  in  false  alarms.  AUG  is  used  to  compare  the  curves  in  Fig.  29 


with  the  test  statistic  maps  displayed  in  Fig.  A  higher  value  indicates  a  better 
quality  test  statistic  map,  meaning  it  is  more  capable  of  providing  higher  accuracy 
change  results.  As  is  seen  in  Table  the  AUG  for  the  nonlinear  optimization  is 
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higher  than  the  uncompensated  CC  and  the  known  mis-registration.  The  fact  that 
the  nonlinear  optimization  approach  provides  better  results  than  the  GLRT  using  the 
known  mis-registration  is  most  likely  due  to  unaccounted  for  mis-registration  artifacts 
in  the  original  data  set. 
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Q2  :  Nonlinear  Optimization  Q1  :  Nonlinear  Optimization 
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Figure  27.  The  four  quadrant  test  statistics  produced  using  the  nonlinear  optimized  test  statistic  for  the  tower  data. 


Minimum  Nonlinear  Optimization 
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Figure  28.  The  uncompensated  CC  test  statistic  (bottom  left),  minimum  nonlinear  optimization  test  statistic  (top),  and  the 
known  mis-registration  (bottom  right)  test  statistic. 


ROC  Curve 


Figure  29.  ROC  curves  for  CC  (solid  line),  proposed  GLRT-based  nonlinear  optimiza¬ 
tion  (dashed  line),  and  the  known  mis-registration  test  statistic  (dotted  line). 


Table  4.  AUC  values  for  the  ROC  curves  displayed  in  Fig. 


29 


Minimum 

Minimum 

CC 

Nonlinear 

Known  Mis- 

Optimization 

registration 

AUC 

0.9973 

0.9980 

0.9974 

4. 5. 1.2  Quadratic  approximation 


The  quadratic  approximation  to  the  quartic  test  statistic  is  compared  to  the  non¬ 
linear  optimization  approach  in  order  to  demonstrate  the  efficacy  of  that  approxi¬ 
mation.  This  is  accomplished  by  investigating  specihc  pixel  classes:  tree,  top-edge 
of  panel,  side-edge  of  panel,  change  pixel,  and  grass.  The  figures  show  the  percent 


error  between  the  test  statistic  {q{x,y,  a,  ^))  from  Eq.  (19)  and  Eq.  (26)  prior  to 
minimization. 
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Fig.  30  shows  the  different  pixel  classes  and  their  spatial  location  in  the  imagery 
used  to  assess  the  efficacy  of  the  approximation.  The  labels  display  the  x—  and  y— 
coordinates,  the  RGB  values,  and  the  color  index  for  each  class  of  pixels.  Fig.  [M] shows 
the  percent  error  between  the  nonlinear  optimization  and  the  quadratic  approximation 
for  Quadrant  II  (where  the  nonlinear  optimization  results  are  regarded  as  ‘truth’).  In 
the  figures,  the  x-axis  and  y-axis  are  a  and  (3  respectively,  while  the  z-axis  corresponds 


to  the  percent  error  between  test  statistics  in  Eq.  (19)  and  Eq.  (26).  Based  on  the 
simulation  results,  tree,  change  pixel,  and  grass  pixel  classes,  it  is  seen  that  the 
minima  of  the  quadratic  approximation  matches  very  closely  to  the  Nelder-Mead 
results  using  the  non-closed  form  solution.  Error  increases  slightly  for  top-edge  and 
significantly  for  side-edge  pixels.  The  results  demonstrate  the  efficacy  of  the  quadratic 
approximation  to  replace  the  nonlinear  optimization,  yielding  a  more  computational 
efficient  option  over  the  Nelder-Mead  numerical  optimization  method.  There  is  a 
30  times  improvement  in  computation  performance  of  the  quadratic  implementation 
over  the  Nelder-Mead  nonlinear  optimization  results. 
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Figure  30.  The  spatial  location  of  investigated  pixels  from  the  image  representing  (top- 
to-bottom)  a  tree,  top-edge  of  panel,  side-edge  of  panel,  change,  and  grass  pixel.  The 
labels  on  the  image  display  the  x—  and  y—  coordinates,  the  RGB  values,  and  the  color 
index  for  each  class  of  pixels. 
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Q2  (72,  109)  :  grass  pixel  02  (94,  26) :  tree  pixel  02  (27. 86)  :  change  pixel 
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nonlinear  optimization  and  the  quadratic  approximation  relative  to  the  nonlinear  optimization  test  statistic  for  Quadrant 


The  test  statistic  for  the  four  quadrants  resulting  from  the  quadratic  approxima¬ 


tion  are  shown  in  Fig.  32  for  =  Uy  =  0.05.  Again,  Quadrant  II  suppresses  the 
mis-registration  effects  observed  on  the  panel  edges  as  seen  by  the  lower  test  statistic. 
Fig,  [^displays  the  uncompensated  CC  result  (bottom  left)  in  comparison  to  the  over¬ 
all  minimum  of  the  quadratic  approximation  (top)  and  the  known  mis-registration 
test  statistic  (bottom  right).  The  true  changes  in  the  scene  are  indicated  with  a  larger 
(white)  test  statistic  while  a  noticeable  smaller  (black)  test  statistic  is  observable  on 
the  panel  edges  while  some  change  errors  occur  from  the  variation  in  shadow  positions 
as  indicated  by  the  shot  noise  appearance  in  the  test  statistic  imagery.  The  ROC  re¬ 
sults  for  the  nonlinear  optimization  and  quadratic  approximation  are  nearly  identical 


(Fig.  34).  This  further  supports  the  use  of  the  more  computationally  efficient  closed 
form  solution  versus  the  nonlinear  form  requiring  numeric  approaches  to  solve.  The 
AUC  value  results  shown  in  Table  highlight  the  similants  of  the  quadratic  approxi¬ 
mation  results  to  the  nonlinear  optimization  results.  Additionally,  both  overcome  the 
known  mis-registration  which  likely  does  not  account  for  mis-registration  artifacts  in 
the  original  data  set. 


Table  5.  AUC  values  for  the  ROC  curves  displayed  in  Fig. 


34 


Minimum 

Minimum 

Minimum 

Nonlinear 

Quadratic 

Known  Mis- 

Optimization 

Approximation 

registration 

AUC 

0.9980 

0.9981 

0.9974 

4.5. 1.3  Mis-registration  estimates 

While  estimating  sub-pixel  shifts  is  not  the  goal  of  this  work,  it  is  instructive  to 
examine  the  statistics  of  a  and  estimates  for  similar  groups  of  pixels  {e.g.,  tree,  top- 
edge  of  panel,  side-edge  of  panel,  change  pixel,  and  grass)  with  both  the  nonlinear 
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GS  ;  Quadratic  Approximation  Q1  :  Quadratic  Approximation 


Figure  32.  The  four  quadrant  test  statistics  produced  using  the  quadratic  approxima¬ 
tion. 

optimization  and  the  quadratic  approximation  for  (t^  =  Cy  =  0.5. 

The  f3  estimate  for  the  top-edge  pixel  is  approaching  the  actual  value  of  0.5  with 
moderate  variance.  The  a  estimate  is  larger,  which  is  expected  since  a  shift  along 
the  edge  has  little  effect.  Based  on  the  known  mis-registration,  (d,  /9)  should  be  (0.5, 
0.5).  Table  provides  the  mean  estimates  of  ot  and  [3  (indicated  with  an  overbar  as 
a  and  /3).  The  d  estimate  for  the  side-edge  pixels  is  near  the  actual  value  of  0.5  with 
moderate  variance.  The  (3  estimate  is  not  close  to  0.5,  presumably  because  a  vertical 
shift  does  not  introduce  a  large  change  for  the  side-edge  pixel.  The  minimum  test 
statistic  for  the  grass  and  tree  pixels  do  not  necessarily  fall  in  the  correct  quadrant. 
This  likely  occurs  because  the  mis-registration-induced  effects  are  small  relative  to 
the  clutter/noise  level  resulting  from  a  lack  of  fine  spatial  content.  Additionally,  the 
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Minimum  Quadratic  Approximation 


.  .1  I 


CC  Known  Mis-registration  Minimum 


Figure  33.  The  uncompensated  CC  (bottom  left),  minimum  quadratic  approximation 
(top),  and  the  known  mis-registration  (bottom  right)  test  statistics  visually  represented. 


clutter/noise  level  impacts  the  a  estimate  on  the  horizontal  edge  and  the  (3  estimate 
on  the  vertical  edge.  The  local  mixing  model  appears  to  weakly  depend  on  the  mis¬ 
registration.  The  algorithm  works  well  with  regard  to  reducing  the  test  statistic  for 
edge  areas  that  are  affected  by  spatial  mis-registration,  even  if  it  does  not  provide  a 
precise  estimate  of  the  mis-registration. 


4.5. 1.4  Covariance  equalization  predictor 


The  experiments  are  repeated  using  the  CE  predictor  to  assess  the  relative  behav¬ 


ior  with  mis-registration  compensation.  As  described  in  Section  2.4.2. 1  CE  estimates 
can  be  substituted  for  CC  estimates.  The  CE  and  CC  test  statistic  images  are  com¬ 
pared  visually  with  and  without  compensation  using  the  quadratic  approximation 
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Figure  34.  ROC  curves  for  test  statistics  based  on  minimum  nonlinear  optimization 
(dashed  line)  and  minimum  quadratic  approximation  (dotted  line). 


- Minimum  Noniinear  Optimization 

- Minimum  Quadratic  Approximation 

Minimum  Known  Mis-registration 


in  Fig.  and  Fig.  The  value  of  the  mis-registration-induced  test  statistic  is 
higher  than  the  true  change,  illustrating  the  impact  mis-registration  has  on  the  un¬ 


compensated  change  detectors  (Fig.  35).  For  both  algorithms,  the  mis-registration 
compensation  signihcantly  suppresses  false  alarms  around  the  panel  edges  (Fig.  [^. 
Comparing  the  corresponding  ROC  curves  in  Fig.  it  is  evident  that  CC  provides 
better  results  than  CE  for  this  case  both  with  and  without  compensation.  Further- 


Table  6.  Mean  estimates  and  variance  of  sub-pixel  shifts  for  Quadrant  II 


Nonlinear  optimization  Quadratic  approximation 


Pixels 

a 

/S 

a 

/S 

Tree 

0.233 

0.145 

0.174 

0.034 

0.303 

0.056 

0.142 

0.055 

Top- Edge 

0.449 

0.139 

0.481 

0.067 

0.531 

0.060 

0.432 

0.028 

Side-Edge 

0.565 

0.138 

0.191 

0.134 

0.661 

0.123 

0.332 

0.044 

Change 

0.426 

0.267 

0.543 

0.167 

0.667 

0.240 

0.833 

0.164 

Grass 

0.302 

0.092 

0.112 

0.014 

0.175 

0.103 

0.109 

0.011 
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more,  the  compensation  algorithm  provides  more  signihcant  improvement  for  CC 
than  CE.  The  AUC  values  are  presented  in  Table  The  AUC  indicates  that  the 
detection  rate  of  the  quadratic  approximation  is  higher  than  the  uncompensated  CC 
or  CE. 


Figure  35.  Visual  results  of  the  uncompensated  CE  change  detector  (left)  and  the  CC 
change  detector  (right). 


Minimum  :  Quadratic  Approximation  (CE)  Minimum  :  Quadratic  Approximation  (CC) 


Figure  36.  The  minimum  quadratic  approximation  results  using  the  CE  estimates  (left) 
and  the  minimum  quadratic  approximation  results  using  the  CC  estimates  (right). 


CE  vs.  CC 


Figure  37.  ROC  curve  for  uncompensated  CE  change  detector  (black  solid)  and  the 
CC  change  detector  (black  dashed)  and  the  minimum  quadratic  approximation  using 
the  CE  (gray  solid)  and  CC  estimates  (gray  dashed). 


Table  7.  AUC  values  for  the  ROC  curves  displayed  in  Fig. 
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Quadratic 

Quadratic 

CC 

approximation 

CC 

CE 

approximation 

CE 

AUC 

0.9972 

0.9981 

0.9909 

0.9954 

4. 5. 1.5  Parameter  sensitivity 


The  GLRT  test  statistic  in  Eq.  (16)  is  based  on  three  model  parameters:  the 


noise  covariance  matrix  C„  and  the  RMS  mis-registration  parameters  and  ay.  The 
value  of  Cn  is  estimated  from  the  difference  image  while  experiment-based  estimates 
for  ax  and  ay  are  required.  To  further  explore  parameter  sensitivity,  the  RMS  mis¬ 
registration  parameter  is  varied  from  0  to  1,  calculating  the  ROC  curve  for  each 
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parameter  value.  The  baseline  CC  algorithm  corresponds  to  =  Cy  =  0.  The 
ROC  curves  in  Fig. [^demonstrate  the  performance  for  the  quadratic  approximation 
algorithm  (with  CC  prediction)  using  Cx  =  <Jy  =  {0.0,  0.05, 0.1,  0.5}.  The  best  results 
are  achieved  with  ax  and  ay  on  the  order  of  0.05  to  0.1.  This  result  indicates  both  the 
importance  of  including  the  prior  distribution  for  the  residual  mis-registration  and  the 
challenge  of  appropriately  setting  the  parameter.  Having  no  prior  distribution  would 
correspond  to  large  ax  and  ay,  which  would  tend  to  poor  performance  as  indicated  in 
(Tx  =  o'y  =  Q  (CC  results)  in  Fig.  38  In  the  next  section,  the  influence  of  parameter 
sensitivity  of  the  image  data  characteristics  and  the  target  mask  employed  in  the 
scoring  is  explored. 


ROC  Curve 


Figure  38.  The  ROC  curve  compares  the  performance  using  the  quadratic  approxi¬ 
mation  for  the  tower  data  with  varied  RMS  mis-registration  parameters  {ox  and  Oy) 
where  CC  {a^  =  ay  =  0)  is  the  solid  line,  ax  =  ay  =  0.05  is  the  dashed  line,  ax  =  ay  =  0.1 
is  the  dotted  line,  and  ax  =  ay  =  0.5  is  the  dash-dot  line. 


4.5. 1.6  Extended  performance  investigation 


In  this  section  the  Mojave  data  set  is  explored  to  provide  additional  insight  on  the 
performance  of  the  proposed  mis-registration  GLRT-based  algorithm.  The  results 
are  examined  against  various  truth  masks  to  demonstrate  the  performance  impact 
of  edge  pixels.  Detection  results  in  the  form  of  ROC  curves  and  AUC  are  used  to 
show  the  performance  of  the  proposed  mis-registration  GLRT-based  algorithm  against 
LCRA.  Upsampling  is  a  pre-processing  requirement  for  sub-pixel  LCRA.  In  order  to 
provide  a  fair  comparison  between  LCRA  and  the  GLRT-based  approach,  images  are 
upsampled  by  two  in  both  the  x  and  y  direction  prior  to  running  the  CD  routine. 

The  proposed  mis-registration  GLRT-based  algorithm  and  LCRA  algorithm  are 
demonstrated  on  both  the  tower  data  and  Mojave  data  sets  for  comparison  purposes. 
Multiple  GLRT  curves  are  generated  by  varying  the  size  of  the  truth  pixel  mask  while 
hxing  the  mis-registration  parameter  ax  =  Cy  =  0.5. 

Fig.  1^  shows  ROC  curves  corresponding  with  an  increasing  erode  of  the  truth 
mask  for  the  tower  data  set.  The  ROC  curves  include  CC,  ax  =  ay  =  0.1,  ax  = 
ay  =  0.5,  and  LCRA.  There  is  one  true  change  in  the  tower  data  set  and  the  ground 
truth  mask  provides  the  spatial  pixel  locations  to  develop  the  ROC  curve  for  every 
threshold.  The  ground  truth  is  defined  at  the  pixel-level  where  the  truth  mask  is 
derived  using  the  CC  test  statistic.  The  target  truth  mask  is  4x4  resulting  in  16 
target  pixels.  A  3  pixel  ignore  mask  extends  the  true  target  mask  for  a  total  size  of 
20x20  pixels.  As  such,  20x20  is  the  maximum  extent  of  the  truth  plus  ignore  pixel 
region.  The  masks  vary  based  on  the  number  of  pixels  associated  with  the  target 
centroid.  The  target  centroid  begins  with  a  4x4  true  change  pixel  mask  and  erodes 
outward  by  1  pixel  in  the  x—  and  y—  direction  toward  the  ignore  mask  edge  per 
ROC  curve  evaluation  iteration.  Thereby,  the  background  pixels  remain  unchanged 
at  66,632  pixels.  The  sizes  of  the  target  centroid  pixels  for  each  ROC  curve  is  4x4, 


6x6,  8x8,  10x10,  and  12x12  (a  pixel  diagram  is  shown  in  Fig.  39  for  reference).  Thus, 
the  ignore  mask  continues  to  shrink  as  the  true  target  centroid  expands.  The  ignore 
mask  begins  with  an  8  pixel  boundary  outward  from  the  4x4  target  pixels  and  shrinks 
by  1  pixel  until  it  reaches  the  12x12  target  with  a  4  ignore  pixel  boundary.  The  ROC 
curves  intersect  thus  the  AUC  is  considered  and  shown  in  Table  Although,  the 
LCRA  tends  to  have  a  higher  AUC  value  for  this  data,  it  must  be  noted  that  the 
AUC  gives  equal  weight  to  the  full  range  of  threshold  values  obscuring  the  fact  that 
the  proposed  GLRT-based  algorithm  performs  better  for  lower  false  alarm  values 
specihcally  those  truth  masks  that  extend  into  the  target  mask  to  reduce  effects  of 
uncertain  edge  pixels. 


Figure  39.  The  diagram  shows  the  target  centroid  surrounded  by  a  pixel  ignore  mask. 
As  the  ROC  performances  are  evaluated  at  each  performance  iteration  the  truth  target 
size  increases  toward  the  edge  of  the  ignore  mask.  The  truth  mask  sizes  of  the  target 
centroid  pixels  for  each  ROC  curve  is  4x4  (red  square),  6x6  (green  square),  8x8  (blue 
square),  10x10  (purple  square),  and  12x12  (orange  square)  for  the  tower  data  set. 
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Table  8.  AUC  values  for  the  ROC  curves  displayed  in  Fig.  40  {tower  data). 


CC 

Quadratic 

approximation 

O' X  eTy  0.1 

Quadratic 

approximation 

Ox  Oy  0.5 

LCRA 

AUC  1 

0.99864 

0.99939 

0.99959 

0.99956 

AUC  2 

0.99875 

0.99949 

0.99955 

0.99974 

AUC  3 

0.99871 

0.99945 

0.99841 

0.99964 

AUC  4 

0.99656 

0.99607 

0.98460 

0.99818 

AUC  5 

0.98135 

0.97968 

0.96089 

0.99137 

Fig.  1^  displays  the  ROC  curves  for  the  Mojave  data  set.  Similar  to  the  tower 
data,  the  ROC  curves  for  CC,  =  o'y  =  0.1,  =  (Ty  =  0.5,  and  LCRA  were 

calculated  for  varying  truth  masks.  The  Mojave  data  contains  two  change  targets 
and  were  again  derived  using  the  CC  test  statistic.  An  ignore  mask  for  both  targets 
extended  3  pixels  past  the  observed  true  target  for  a  total  size  of  24x24  pixels.  Both 
target  truth  masks  varied  at  the  same  rate  per  evaluation.  The  target  pixels  begin 
with  a  2x2  square  true  change  region  before  expanding  outward  toward  the  edge  of 
the  ignore  mask.  The  background  pixels  account  for  a  total  of  196,860  pixels.  The 
sizes  of  the  target  centroid  pixels  for  each  ROC  curve  calculation,  starting  with  the 
upper  left  plot  of  Fig.  is  2x2,  6x6,  10x10,  14x14,  and  18x18.  The  ignore  mask 
starts  with  an  11  pixel  boundary  and  reduces  to  a  3  pixel  boundary  for  the  hnal  ROC 

are  provided  in  Table 


curve.  The  AUC  values  for  the  curves  in  Fig.  41 
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Table  9.  AUC  values  for  the  ROC  curves  displayed  in  Fig.  41  {Mojave  data). 


CC 

Quadratic 

approximation 

ax  ay  0.1 

Quadratic 

approximation 

ax  0.5 

LCRA 

AUC  1 

0.99256 

0.99666 

0.99800 

0.99793 

AUC  2 

0.99240 

0.99679 

0.99816 

0.99824 

AUC  3 

0.99307 

0.99735 

0.99862 

0.99863 

AUC  4 

0.99362 

0.99789 

0.99916 

0.99913 

AUC  5 

0.99279 

0.99512 

0.97564 

0.95889 

After  observing  the  {ax,  ay)  parameter  sensitivity  with  the  varying  truth  analysis, 
the  change  detection  performance  is  found  to  vary  somewhat  based  on  user-defined 
parameters  for  the  underlying  prior  probability  distribution  for  the  residual  misreg¬ 
istration.  While  this  supports  the  importance  of  including  a  prior  distribution  in 
the  model,  it  motivates  further  work  to  better  understand  this  parameter  sensitivity. 
This  specihc  experiment  demonstrates  that  the  algorithmic  performance  is  largely 
based  on  coverage  of  the  truth  mask.  The  conclusion  of  which  algorithm  has  superior 
performance  can  change  depending  on  the  truth  mask  used  in  the  analysis.  For  the 
middle  truth  mask,  the  proposed  mis-registration  GLRT-based  approach  and  LCRA 
are  comparable. 
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LCRA  (dashed  line). 


The  (a^,  ay)  parameter  tuning  results  for  the  propsed  GLRT-based  algorithm  on 
the  interpolated  data  justify  further  exploration  of  the  parameter  and  its  impact 
on  performance.  The  RMS  mis-registration  parameter  is  varied  from  0  to  1  for  the 
Mojave  data  set,  calculating  the  ROC  curve  for  each  parameter  value.  The  truth  mask 


in  Fig.  42  (left)  contains  one  change  target  as  dehned  at  the  pixel-level  and  derived 
using  the  CC  test  statistic.  There  is  a  border  surrounding  the  truth  target  that  is 
ignored  during  evaluation  to  eliminate  uncertainty  along  the  target  edges.  The  white 
target  pixels  are  a  6x6  square  consisting  of  36  pixels  and  the  boundry  pixels  radiate 
outward  7  pixels  in  each  direction  creating  a  border  of  364  pixels  (or  a  20x20  square 
region  excluding  the  target  pixels).  The  total  number  of  background  pixels  is  66,632. 


The  baseline  CC  algorithm  corresponds  to  a^  =  ay  =  0.  The  ROC  curves  in  Fig.  42 


(right)  demonstrate  the  performance  for  the  quadratic  approximation  algorithm  (with 
CC  prediction)  using  a^  =  ay  =  {0.0,0.05,0.1,0.5}  on  the  tower  data.  The  best 
results  on  the  upsampled  imagery  are  achieved  with  a^  and  ay  are  on  the  order  of 
0.5  and  0.1.  In  the  non  upsampled  domain,  the  results  correspond  to  0.25  and  0.05 

even 


respectively.  The  results  are  consistent  with  results  presented  in  Section  4. 5. 1.5 


though  the  truth  mask  may  be  somewhat  different.  The  variation  in  results  is  due  to 
a  difference  in  edge  pixels  between  the  upsampled  and  original  imagery. 

In  addition,  the  Mojave  data  set  was  used  to  explore  the  variation  of  performance 


with  the  mis-registration  parameters.  The  truth  mask  in  Fig.  43  (left)  contains  two 
change  targets,  both  have  ignore  pixels  surrounding  the  target  edges.  The  left  and 
right  changes  include  an  8x8  target  square  (64  pixels)  with  8  ignore  boundry  pixels 
(24x24=576  square  region  excluding  the  target  pixels).  The  total  number  of  back¬ 
ground  pixels  is  196,860.  The  ROC  curves  are  generated  by  varying  the  (a^,  ay) 


parameter  and  are  shown  in  Fig.  43  (right).  The  best  results  are  achieved  with  the 
upsampled  Mojave  data  set  are  on  the  order  of  0.25  and  0.05  in  the  original  non 
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upsampled  space  (Section  4. 5. 1.5).  The  conclusion  is  the  best  performing  parameter 
range  is  consistent  through  repeated  results  on  different  imagery  and  varying  truth 
masks. 


Truth  map 


ROC  Curve 


Figure  42.  Parameter  investigation  for  {ax,<Ty)  with  the  tower  data.  The  truth  mask 
(left)  displays  target  pixels  in  white  and  boundry  pixels  in  gray  which  are  ignored  during 
evaluation.  The  ROC  curve  (right)  compares  the  performance  using  the  quadratic 
approximation  with  varied  RMS  mis-registration  parameters  {ax,a-y)  where  CC  {(Jx  = 
ay  =  0)  is  the  solid  line,  ax  =  ay  =  0.05  is  the  dashed  line,  ax  =  ay  =  0.1  is  the  dotted  line, 
and  ax  =  ay  =  0.5  is  the  dash-dot  line. 


4.5.2  Parallax-compensation  algorithm 

This  section  discusses  the  results  for  the  parallax  compensation  algorithm  derived 
in  Section  |3.2[  The  primary  experimental  results  are  promising  by  visual  inspection 
as  any  potential  quantitative  assessment  would  be  difficult  to  derive  with  lack  of  truth 
data.  For  the  synthetic  DIRSIG  data  set,  the  assessment  is  performed  visually  and 
good  change  detection  is  obtained  while  minimizing  parallax  errors. 

4. 5. 2.1  Synthetic  DIRSIG  data 

This  section  shows  that  the  correspondence  matching  pixel  shift  estimate  from 
Section  |3.2|  provides  a  good  approximation  of  the  parallax  direction  (denoted  here  as 
the  computed  parallax  direction).  The  results  for  the  manually  estimated  parallax 
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Truth  map 


ROC  Curve 


Figure  43.  Parameter  investigation  for  {ax,o'y)  with  the  Mojave  data.  The  truth  mask 
(right)  displays  target  pixels  in  white  and  boundry  pixels  in  gray  which  are  ignored 
during  evaluation.  The  ROC  curve  (left)  compares  the  performance  using  the  quadratic 
approximation  with  varied  RMS  mis-registration  parameters  {ax,a-y)  where  CC  {(Jx  = 
ay  =  0)  is  the  solid  line,  ax  =  ay  =  0.05  is  the  dashed  line,  ax  =  ay  =  0.1  is  the  dotted  line, 
and  ax  =  ay  =  0.5  is  the  dash-dot  line. 


direction  (90°  as  computed  directly  from  the  imagery)  are  shown  visually  in  Fig.  44 


The  parallax  error  parameter  (cTp)  is  varied  while  the  mis-registration  parameter  is 
held  constant  (cxm  =  0.5)  while  searching  over  (m,h)  =  (— 1,  {— 8,  — 7...7,  8})  and 
{m,n)  =  (0,  {—8,  —7.. .7,  8})  regions  of  whole  pixel  shifts  (double  in  extent  in  the  n 

is 


direction  relative  to  the  known  maximum  parallax  shift).  The  hrst  entry  in  Fig.  44 


the  baseline  CC  and  the  visual  result  for  increasing  Up  follow.  According  to  Fig.  44 


as  cTp  progresses  from  low  (too  severe  a  penalty  function)  to  high  (too  lax  a  penalty 
function),  best  performance  occurs  at  cTp  =  4.0  pixels  (approximately  the  true  paral¬ 
lax). 


Fig.  45  shows  the  per  pixel  mapping  of  the  minimum  value  of  g(  ■ )  in  Eq.  (41) 
resulting  from  the  {rh,  h)  search  window  between  a  hxed  pixel  in  the  reference  image 
to  the  search  area  in  the  test  image  for  a  hxed  value  of  cXp.  For  the  ground  plane  where 
no  parallax  exists,  values  of  {fh,h)  should  be  approximately  (0,0).  This  observation 
is  shown  in  the  visual  (m,  n)  maps  as  the  ground  plane  values  are  represented  by 
n  =  0.  For  buildings  that  exhibit  parallax,  values  of  {rh,n)  should  approximate  the 
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parallax  shift  (in  pixels)  along  edges  in  the  parallax  direction.  The  n  values  should 


correspond  roughly  to  the  truth  map  shown  in  Fig.  23  The  value  of  Cp  stabilizes  at 
approximately  cTp  =  4  as  shown  in  the  (m,  h)  mapping  plots  where  the  building  edges 
take  on  a  vertical  value  shift  of  h  according  to  the  truth  map. 


Fig.  46  displays  the  results  using  the  computed  parallax  direction  attained  by 


Eq.  (46)  versus  the  manually  estimated  parallax  direction.  The  computed  parallax 
angle  is  0  =  86.6°  which  is  very  close  to  the  manually  estimated  parallax  direction 


of  90°.  The  minimum  selection  map  (Fig.  47)  is  now  well-structured  and  the  h  =  —4 
pixels  on  the  building  edges  are  easily  identihed.  Results  between  the  manually 
estimated  0  and  the  computed  value  0  show  signihcant  differences  in  structure.  This 
indicates  some  sensitivity  in  the  value  of  0  during  the  optimization.  The  difference 
between  the  manual  estimate  and  the  computed  are  likely  due  to  a  slight  sub-pixel 
parallax  shift  attained  during  the  registration  process  of  the  temporal  images.  The 
results  show  there  is  little  sensitivity  to  the  parallax  parameter  (up)  and  for  values  of 
CTp  greater  than  the  manually  estimated  parallax  direction,  there  is  no  penalty  paid. 
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little  sensitivity  to  the  parallax  parameter,  ap  and  at  a  ap  greater  than  the  parallax,  there  is  no  penalty  paid  for  improvement. 


The  calculated  synthetic  ground  truth  mask  described  in  Section  |4.2.2.1  provides 
further  insight  into  the  performance  of  the  algorithm.  The  truth  parallax  class  pixels 
mean  and  variance  were  calculated  as  1.9124  and  0.7228,  respectively,  as  dehned  in 


the  truth  mask  in  Fig.  23  The  computed  parallax  direction  maps  computed  via  the 
correspondence  algorithm  shown  in  Fig.  |47|  were  limited  to  n  as  the  parallax  shift  is 
vertical  in  order  to  compare  the  parallax  shift  in  pixels  to  the  truth  parallax  shift.  The 


visual  results  are  displayed  in  Fig.  The  hrst  image  corresponds  to  the  calculated 
parallax  truth  based  on  the  heights  of  the  buildings.  The  middle  of  the  buildings 
are  then  catergorized  as  ignore  pixels  unless  the  pixels  are  along  the  edges  signifying 
parallax  pixels.  The  subsequent  images  are  the  visual  results  cooresponding  to  the  n 
shift  as  selected  by  the  proposed  parallax-compensation  algorithm  as  (jp  is  increased. 
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Figure  48.  The  first  image  corresponds  to  the  calculated  parallax  truth  based  on  the  heights  of  the  buildings.  The  middle 
of  the  buildings  are  then  catergorized  as  ignore  pixels  unless  the  pixels  are  along  the  edges  signifying  parallax  pixels.  The 
subsequential  images  are  the  results  cooresponding  to  the  h  shift  as  selected  by  the  proposed  parallax-compensation  algorithm 
as  (T„  is  increased. 


Table  10  tabulates  the  mean  and  variance  for  each  Up  ranging  from  zero  to  one  (for 
am  =  0.5)  for  the  ground  plane  pixels  and  the  parallax  pixel  results  in  the  minimum 


selection  map  (Fig.  48).  Accordingly,  the  ground  plane  pixel  statistics  are  close  to 
zero  which  is  anticipated  but  as  Up  increases  so  does  the  error.  The  parallax  mean 
statistics  get  closer  to  the  actual  known  value  as  Up  increases.  The  variance  of  the 
parallax  pixels  also  increases  which  increases  the  uncertainty  of  the  mean  estimates. 
In  addition,  the  ground  plane  mean  and  variance  diverge  from  the  calculated  value 
becoming  less  accurate  while  Up  increases  beyond  the  true  parallax  shift.  Both  cases 
result  in  the  increased  clutter  appearance  in  Fig.  The  impact  is  that  true  change 
pixels  could  be  lost  as  cXp  grows  too  large. 


105 


Table  10.  Statistics  of  the  detection  statistic  q{  • )  for  the  parallax  and  ground  plane  pixels. 


(Tp 

•) 

d( 

■) 

0 

ground 

parallax 

0.504565 

0.880782 

ground 

parallax 

0.292905 

0.704813 

1 

ground 

parallax 

0.753677 

1.405863 

ground 

parallax 

0.555484 

1.276498 

2 

ground 

parallax 

1.043424 

1.629967 

ground 

parallax 

0.81114 

1.510966 

3 

ground 

parallax 

1.155527 

1.686645 

ground 

parallax 

0.888017 

1.555591 

4 

ground 

parallax 

1.203219 

1.715961 

ground 

parallax 

0.917929 

1.593324 

5 

ground 

parallax 

1.227652 

1.725081 

ground 

parallax 

0.935127 

1.598425 

6 

ground 

parallax 

1.242345 

1.72899 

ground 

parallax 

0.944134 

1.601864 

7 

ground 

parallax 

1.251017 

1.731596 

ground 

parallax 

0.949046 

1.601967 

8 

ground 

parallax 

1.256184 

1.736156 

ground 

parallax 

0.951562 

1.607655 

9 

ground 

parallax 

1.259919 

1.738111 

ground 

parallax 

0.953887 

1.614551 

10 

ground 

parallax 

1.262952 

1.741368 

ground 

parallax 

0.954663 

1.618203 

Further,  performance  is  visualized  through  a  sequence  of  magnihed  test  statistic 
images  (here,  the  magnihcation  refers  to  zooming  in  on  the  building  of  interest  vice 
processing  on  the  test  statistic  itself)  of  a  few  of  the  buildings  from  the  CC  CD  with 


Up  =  0  to  Up  =  10.  For  reference.  Fig.  49  shows  the  buildings  in  the  scene  designated 


with  a  numeric  label.  The  results  for  building  label  3  are  shown  in  Fig.  where  it 
is  evident  that  the  increase  in  Up  decreases  the  test  statistic  surrounding  the  build¬ 
ing  which  is  caused  by  parallax.  These  magnihed  test  statistic  images  are  for  the 
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computed  parallax  direction  derived  via  correspondence  matching  (0  =  86.6°).  The 


corresponding  histograms  are  shown  in  Fig.  which  shows  the  distribution  changes 
for  the  parallax  pixels  as  cXp  approaches  the  manually  estimated  parallax  direction. 
The  histograms  show  increasing  pixel  occurrences  for  the  lower  brightness  values  (he., 
pixel  has  been  re-classi£ed  to  a  lower  brightness  value)  as  Up  steps  to  the  next  incre¬ 
ment.  In  other  words,  the  histograms  represent  the  false  alarm  mitigation  activity, 
a  larger  test  statistic  pixel  value  is  mitigated  to  a  smaller  test  statistic  pixel  value 
as  (Tp  increases  toward  the  manually  estimated  parallax  direction.  This  is  repeated 
for  building  labeled  6  to  emphasize  the  reduction  in  parallax-induced  false  alarms 


(Fig.  52  and  Fig.  53). 


Figure  49.  RGB  of  DIRSIG  image,  where  buildings  are  labeled  1  —  7. 
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Minimum  (o  =  0)  Minimum  (a  =  0.5)  Minimum  (a  =  1) 
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Figure  50.  Magnified  test  statistic  results  of  building  3  to  visualize  reduction  in  false  alarms  due  to  parallax. 


Minimum  (cr  =  0)  Minimum  (a  =  0.5)  Minimum  (ci  =  1) 
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Figure  51.  Test  statistic  histograms  for  building  3  for  increasing  Up.  The  histograms  show  a  redistribution  of  parallax  pixels  in 
Fig.  50  as  ffp  increases. 
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Figure  52.  Magnified  test  statistic  results  of  building  6  to  visualize  reduction  in  false  alarms  due  to  parallax. 
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Figure  53.  Test  statistic  histograms  for  building  6  for  increasing  Up.  The  histograms  show  a  redistribution  of  parallax  pixels  in 
Fig.  52  as  ap  increases. 


Although  ROC  curve  analysis  is  preferred,  performance  would  be  hampered  by 
other  spurious  false  alarm  sources.  Additional  false  alarms  sources  not  caused  by 
parallax  can  affect  the  reliability  of  the  change  detector  system  and  are  in  influenced 
by  factors  such  as  shadows,  illumination  differences,  and  atmospheric  effects. 

In  addition  to  the  visual  magnified  test  statistic  images,  the  SCR  is  used.  The 
SCR  is  computed  individually  for  each  building  defined  as 


SCR  = 


/^Change  target  /^Building 
^Building 


(60) 


where  /ichange  target  snd  /iBuilding  Siie  the  mean  of  the  change  target  and  the  mean  of  the 
building  of  interest,  respectively,  and  aBunding  is  the  standard  deviation  of  the  building 


of  interest.  The  plot  shown  in  Fig.  54  is  the  SCR  versus  ap  and  demonstrates  how  the 
clutter  from  the  edge  artifacts  is  reduced  with  an  increase  of  ap  resulting  in  a  higher 
signal.  The  SCR  plot  suggests  for  a  large  ap  there  is  no  penalty  paid  for  this  parallax 
compensation  improvement.  This  is  likely  because  the  change  target  is  spectrally 
different  than  the  local  surroundings.  Otherwise,  the  SCR  would  start  to  decrease  at 
some  point  or  perhaps  not  increase  as  ap  increases  because  the  algorithm  would  find 
a  local  match  for  the  true  change  target.  Contrary  to  the  SCR  plot,  the  statistics 
(Table  [Io|)  and  the  visual  interpretation  images  indicate  that  a  ap  value  greater  than 
the  true  parallax  would  result  in  missed  true  change  pixels. 


4.6  Summary 

Incorporating  a  spatial  mis-registration  model  into  a  GLRT-based  CD  methodol¬ 
ogy  leads  to  a  fourth  order  polynomial  test  statistic  that  must  be  numerically  mini¬ 
mized  to  detect  changes  in  hyperspectral  imagery.  Results  using  controlled  imagery 
with  synthetically  introduced  mis-registration  show  that  the  approach  is  able  to  sup- 
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- Building  1 

Building  2 

- Building  3 

- Building  4 

Building  5 
Building  6 
- Building  7 


Figure  54.  SCR  plot  for  Buildings  1-7  as  referenced  in  Fig.  49  of  the  synthetic  DIRSIG 


data.  This  plot  shows  there  is  no  penalty  for  the  parall2Lx  compensation  improvement 
as  a„  increases  because  the  SCR  does  not  decrease. 


press  mis-registration-introduced  false  alarms  due  to  high  spatial  frequency  content 
as  seen  at  edges  caused  by  larger  contrast  differences  {e.g.,  edges  of  the  panel  and 
hne  spatial  structure  with  tree  leaves).  While  the  nonlinear  optimization  is  numer¬ 
ically  complex,  a  closed-form  quadratic  approximation  provides  nearly  identical  CD 
performance  while  running  30  times  faster.  While  the  mis-registration  compensation 
approach  provides  a  beneht  for  both  CC  and  CE  change  predictors,  a  much  more 
substantial  improvement  in  the  presence  of  residual  mis-registration  is  achieved  using 
the  CC  method. 

The  CD  performance  is  found  to  vary  somewhat  based  on  user-dehned  parameters 
for  the  underlying  prior  probability  distribution  for  the  residual  mis-registration.  The 
best  performance  appears  to  be  data-dependent.  While  this  supports  the  importance 
of  including  a  prior  distribution  in  the  model,  a  unique  contribution  of  our  approach 
over  prior  methods,  it  also  motivates  further  work  to  better  understand  this  parameter 
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sensitivity. 

The  largest  contribution  of  our  generalized  CD  approach  is  that  it  can  be  ex¬ 
tended  to  incorporate  more  complicated  prior  probability  distributions  for  residual 
mis-registration  and  serves  as  a  starting  point  to  account  for  other  known  image 
acquisition  and  preprocessing  artifacts.  In  particular,  current  hyperspectral  CD  al¬ 
gorithms  (CC  and  CE)  are  leveraged  in  conjunction  with  bilinear  interpolation  and 
stereo  correspondence  to  develop  an  algorithm  for  performing  hyperspectral  CD  to 
address  sub-pixel  mis-registration  and  parallax  effects.  The  experiments  clearly  show 
that  parallax  mitigation  using  correspondence  matching  gives  superior  performance 
when  compared  to  the  manual  parallax  estimate  approach.  The  visual  detections  ex¬ 
emplify  the  ability  of  the  hyperspectral  algorithm  to  reduce  parallax  by  eliminating 
the  strong  test  statistic  at  building  or  structure  edges.  In  addition  to  the  visual  perfor¬ 
mance  of  the  algorithm,  the  test  statistic  histograms  and  SCR  computations  provide 
evidence  of  the  reduction  in  parallax-induced  edge  artifacts  in  HSCD  performance. 
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V.  Conclusion 


In  this  final  chapter,  a  summary  of  the  research  is  provided.  Additionally  research 
contributions  are  highlighted  and  the  possibility  of  future  work  discussed. 

5.1  Summary 

A  summary  of  change  detection  (CD)  processes  was  discussed  and  a  framework  was 
outlined  to  understand  the  design  space  of  CD  applications  in  Chapter  2.  A  detailed 
mathematical  formulation  of  the  algorithms  developed  was  presented  in  Chapter  3. 
The  proposed  mis-registration  algorithm  compensates  for  mis-registration  in  the  CD 
process  based  on  minimizing  the  detection  test  statistic  derived  using  the  generalized 
likelihood  ratio  test  (GLRT).  It  was  demonstrated  that  the  GLRT  formulation  results 
in  a  closed  form  solution  for  the  test  statistic  minimized  relative  to  the  unknown  local 
shift  parameters.  The  mathematical  formulation  led  to  a  fourth-order  polynomial 
solution  that  was  minimized  numerically.  By  performing  a  second  order  Taylor  series 
expansion  of  the  test  statistic,  a  quadratic  approximation  was  developed  leading  to  a 
less  computationally-complex  closed-form  analytical  solution.  Moreover,  the  GLRT 
formulation  provided  a  gateway  to  address  more  complex  false-alarm  sources  typically 
found  in  persistent  surveillance  scenarios.  To  that  end,  the  mis-registration  GLRT 
test  statistic  was  extended  to  address  parallax  compensation  by  rotating  the  axis  along 
the  parallax  direction  to  account  for  the  relative  orientation  between  two  consecutive 
image  views. 

In  Chapter  4,  experimental  results  showed  that  the  methods  developed  in  this 
research  achieve  good  performance  when  applied  to  real  or  synthetic  hyperspectral 
data  sets.  This  was  demonstrated  by  comparing  performances  with  those  reached  by 
two  standard  anomalous  CD  methods,  chronochrome  (CC)  and  covariance  equaliza- 
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tion  (CE).  In  the  case  of  the  mis-registration  GLRT,  performance  was  additionally 
compared  with  the  local  co-registration  adjustment  (LCRA)  algorithm.  The  results 
indicated  far  superior  performance  relative  to  CC  and  CE,  and  comparable  perfor¬ 
mance  relative  to  LCRA.  An  important  beneht  to  the  approach  developed  in  this 
dissertation,  however,  is  the  ability  to  extend  it  to  other  more  complex  signal  pro¬ 
cessing  functions,  an  attribute  not  exhibited  by  LCRA. 

Change  detector  performance  was  assessed  through  a  visual  analysis  and  quanti¬ 
tatively  using  the  receiver  operating  characteristic  curves  (ROC)  and  the  area  under 
the  ROC  curve  (AUC)  subject  to  the  availability  of  detailed  truth  information.  Ex¬ 
periments  showed  that  signihcant  performance  improvement  is  accomplished  in  the 
presence  of  mis-registration  for  CC  (AUCqlrtcc  ~  0.9981  and  AUCcc  =  0.9972)  and 
CE  (AUCqlrtce  =  0.9954  and  AUCce  =  0.9909).  The  CLRT  showed  similar  per¬ 
formance  over  LCRA  for  sub-pixel  shifts.  Experiments  clearly  showed  that  parallax 
artifacts  in  the  detection  statistic  are  reduced  compared  to  CC.  LCRA  is  not  designed 
to  address  this  type  of  structured  mis-registration  and  was  therefore  not  evaluated 
against  it. 

5.2  Research  contributions 

The  contributions  of  this  work  included  practical  and  theoretical  aspects.  From 
the  practical  point  of  view,  the  CD  framework  helps  the  CD  researcher  identify  re¬ 
search  areas  of  interest  and  helps  the  analyst  build  functional  CD  systems.  From  the 
theoretical  point  of  view,  the  mis-registration  CLRT  test  statistic  provided  a  gateway 
to  more  complex  signal  processing  tasks  related  to  CD  and  the  quadratic  approxima¬ 
tion  to  the  CLRT  test  statistic  provided  a  practical  implementation  of  that  approach. 
Extending  the  mis-registration  CLRT  test  statistic  to  include  parallax  compensation 
for  hyperspectral  change  detection  provided  a  novel  contribution  in  that  such  an 
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algorithm  has  yet  to  be  published  in  the  literature. 


5.3  Future  work 

Although  the  research  presented  in  this  dissertation  document  is  highly  successful 
and  presents  novel  approaches  to  reducing  complex  false- alarm  sources  in  CD,  there 
are  possible  areas  for  future  work.  These  are  outlined  briefly  in  this  section. 

As  shown  in  Chapter  4,  excellent  false-alarm  mitigation  results  are  achieved  us¬ 
ing  (SIFT)  features  for  the  correspondence  matching  in  the  parallax  compensation 
algorithm.  However,  other  correspondence  algorithms  exist  and  may  improve  perfor¬ 
mance  by  providing  more  reliable  feature  matching  between  the  test  and  reference 
images.  Possible  correspondence  approaches  include  contrast  invariant  matching  1^2] 
and  iterative  approaches  [93].  In  doing  so,  however,  the  correspondence  matching 
approach  may  be  more  complex.  Emerging  new  research  areas  such  as  computational 
photography  give  rise  to  novel  applications  for  correspondence  problems  that  may  be 
applicable  to  the  current  body  of  research  [M]. 

A  critical  component  in  the  calculation  of  the  parallax  direction  is  the  accuracy 
in  mapping  corresponding  pixels  between  image  pairs.  Higher  accuracy  of  the  cor¬ 
respondences  results  in  a  more  accurate  estimate  of  the  parallax  direction.  Another 
consideration  in  an  urban  surveillance  scenario  is  that  building  occlusions  hinder 
the  matching  of  potentially  good  image  features.  Investigation  of  mismatches  or  ap¬ 
proaches  accounting  for  occlusion  pixels  may  uncover  a  more  accurate  correspondence 
matching  method.  It  is  recognized  in  the  literature  that  multiple  sources  of  correspon¬ 
dence  matching  errors  exist  [9l] .  These  error  sources  include  left-right  consistency  and 
the  flatness  of  the  cost  curve  [2S|.  These  errors  are  mitigated  by  using  the  matching 
cost  with  heuristics  rules  [95] .  Occlusion  is  a  relevant  problem  to  mismatching  errors 
and  there  exist  three  classes  of  algorithms  for  handling  occlusion;  methods  that  de- 
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tect  occlusion,  methods  that  reduce  sensitivity  to  occlusion,  and  methods  that  model 
the  occlusion  geometry  [96].  Correspondence  robustness  and  optimization  matching 
should  be  considered  in  future  work.  The  matching  algorithm  may  favor  a  target 
data  set,  for  example,  hyperspectral  data,  multi-spectral  data,  or  intensity  (gray¬ 
scale)  data.  Furthermore,  the  user  must  decide  whether  or  not  any  accuracy  gains  in 
more  complex  correspondence  matching  approaches  are  worth  the  potential  increase 
in  computational  time. 

The  parallax-compensation  algorithm  requires  a  (m,  n)  search  space  dehned  by 
the  user.  Future  research  can  focus  on  automatizing  the  (m,  n)  search  space  through 
the  estimated  parallax  direction.  For  example,  targeting  a  larger  search  along  the  (j) 
direction  and  limiting  the  orthogonal  direction  to  a  smaller  search  for  efficiency. 

The  CD  performance  is  found  to  vary  somewhat  based  on  user-defined  parame¬ 
ters  for  the  underlying  prior  probability  distribution  for  the  residual  mis-registration. 
In  fact,  the  best  performing  parameter  range,  as  determined  through  experimenta¬ 
tion,  was  found  to  be  consistently  lower  than  the  known  mis-registration.  While  this 
supports  the  importance  of  including  a  prior  distribution  in  the  model  (a  unique 
contribution  to  this  work)  it  also  motivates  further  work  to  better  understand  the 
relationship  of  the  model  parameters  to  the  known  mis-registration.  Finally,  the  ad¬ 
dition  of  more  extensive,  real,  truthed,  hyperspectral  data  to  exercise  and  test  would 
refine  the  approach.  This  requires  such  data  to  be  available,  which  is  currently  a 
major  shortcoming. 
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